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ABSTRACT 

An automatic control system capable of controlling an 
unknown nonlinear system is designed using a direct adaptive 
control scheme, implemented with a Hopfield network. The 
application of this method to an arbitrary system is discussed 
in detail and three specific simulation studies are included. 
These studies include the implementation of the Hopfield 
network based direct adaptive control system to a linear 
system, an inverted pendulum, and a nonlinear model of the NPS 
Autonomous Underwater Vehicle (AUV) with six degrees of 
freedom. The AUV simulation includes a three dimensional 
trajectory following algorithm and shows the ability of the 
Hopfield network to adapt to simultaneous ordered changes in 
the AUV's depth, speed, and course. 

Additionally, an analog circuit design is proposed which 
implements the automatic control scheme without the support of 
a microprocessor. The circuit was set up in SPICE and the 
simulation results as well as some limitations of the analog 
circuit implementation of the Hopfield network are presented. 
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I. INTRODUCTION 



Apart from a few particular cases, no general theory 
exists for the control of nonlinear systems. The simplest 
method for the development of a control system for a nonlinear 
plant is to linearize the nonlinear plant around an operating 
point and derive a linear controller with classical control 
methods. These static control algorithms can control 
nonlinear systems at specific operating points, but may be 
unstable at others. Simple nonlinear systems such as the 
inverted pendulum may have stable linearizations at some 
operating points and be unstable at others. These facts make 
the control of nonlinear systems difficult. 

Adaptive control algorithms hold promise in the control of 
nonlinear systems. Because an adaptive control system changes 
in response to changes in the system, it is able to control 
many nonlinear plants. This thesis investigates the use of a 
Hopfield net for direct adaptive control of a nonlinear 
system. 

An interesting application is the adaptive control of an 
Autonomous Underwater Vehicle (AUV) . An AUV is an unmanned 
submersible vehicle designed to operate independently of human 
interaction or support. As such, it must be capable of 
responding to changing, dangerous, or unpredictable conditions 
much as a manned underwater vehicle would. While many of the 
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AUV’s higher level functions, such as path planning and 
obstacle avoidance, use artificial intelligence techniques to 
cope with different scenarios, the lower level control remains 
unsolved. 

The dynamics of the AUV are highly non-linear and are not 
easily rendered into a satisfactory linear form. These non- 
linearities are particularly evident as the vehicle changes 
speeds. Schwartz investigated the use of recursive least 
squares (RLS) and an adaptive pole placement scheme to control 
the AUV at a given speed [Ref. l:p. 63]. Although this method 
produced functional results, the intensive calculations 
required by the RLS algorithm would further burden the already 
heavily loaded microprocessor. Since robustness to changing 
plant parameters (due to changing environmental conditions or 
damage) is required, an adaptive controller must be used. As 
the processor aboard is already occupied by the path planning, 
sensor data processing, and navigation software, a solution to 
the nonlinear automatic control problem was sought that 
produces a controller that demands less processor time. 

Neural networks offer a potential solution to this problem 
as there is great promise in the implementation of neural nets 
in analog hardware. The Hopfield network in particular is 
easily realized in analog circuitry. An adaptive controller 
designed using a Hopfield network realized in analog hardware 
would not overload the on-board processor and yet would 
provide the necessary robustness for the control of the AUV. 
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The goal of this thesis is to develop an adaptive control 
algorithm based on the Hopfield network and to propose a 
design in analog hardware. 

This thesis is organized in five sections. Chapter II 
describes the direct adaptive control algorithm used 
throughout this thesis. Chapter III contains a description of 
the Hopfield network and its application to direct adaptive 
control. Chapter IV consists of three studies of the 
implementation of the Hopfield network based direct adaptive 
control scheme. Chapter V describes a possible analog circuit 
implementation of this control scheme and discusses some of 
the problems associated with it. Chapter VI provides a summary 
of the results of this work and points out several areas for 
further study. 
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II. ADAPTIVE CONTROL 



Adaptive control is a method by which a controller is 
adapted to control an uncertain system in a dynamic operating 
environment. Two major classes of adaptive control exist: 
direct and indirect adaptive control. Direct adaptive control 
is characterized by the direct determination of the control 
parameters from input and output data collected from the 
system. Indirect adaptive control is a two stage process. 
First, system identification techniques are used to obtain a 
model of the system and then standard control techniques are 
used to calculate a controller for the estimated model [Ref. 
2:p. 48]. Indirect adaptive control is generally slower than 
direct adaptive control and requires greater hardware support 
and/or computational effort. In this thesis only direct 
adaptive control is considered. 

A. PARTIAL STATE REPRESENTATION 

In order to proceed with the derivation of the direct 
adaptive control algorithm, it is necessary to introduce an 
alternate representation of a system called the partial state 
representation [Ref 3:p. 209]. The system 

y(t) = £lfiu(t) (2-1) 
A(s) 



4 



where s can be interpreted either as the differential operator 
or the complex variable of the Laplace transform and 




can be broken into two components as shown in Figure 2-1. 



The intermediate state z(t ) is called the partial state and 
the system of equation (2-1) is equivalent to 



The partial state representation is useful in the derivation 
of the direct adaptive control algorithm. 

B. DIRECT ADAPTIVE CONTROL 

As stated earlier, direct adaptive control uses the input 
signal to a system and the output signal from a system to 
directly determine suitable control parameters. For this work 
a pole placement scheme is employed, meaning that the system 




y(t) 



Figure 2-1 - Partial State Representation 



A[s) z (t) = u(t) 

y(t) =B(s)z(t) 



(2-3) 
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output y ( t ) is 
controlled to react to 
the reference signal 
v(t) as does the 
reference model with the 
transfer function 
l/p'(s). A block diagram 
of a direct adaptive 
control system is shown 
in Figure 2-2. In this diagram, the control parameter 
identifier receives data from the input and output of the 
system and uses this data to modify the controller. 

Assuming the unknown plant can be modeled as a piecewise 
linear system, then it may be described at any given operating 
point by the linear differential equation: 

y(t) = (2-4) 

P(s) 

where p(s) is an N * h order monic polynomial and r(s) is an M th 
order stable polynomial, meaning that the roots of r(s) are 
all in the left half-plane. Rewriting equation (2-4) in 




Figure 2-2 - Direct Adaptive 

Control System 
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partial state form yields: 

p(s)z(t) = u(t) 2 _ 5 

y(t) = r (s)z(t). 

The goal of the controller is to track the output of a 
reference model driven by an external input v(t) , specifically 

yjt) = — - — v(t) (2-6) 

P’(s) 

where p'(s) is the characteristic polynomial of the reference 
model and y m (t) is the reference model output. From the pole 
placement problem the control input has the following 
structure 

u(t) * ^ifly(t)+^lflu(t)+g -v(t) (2-7) 

q(s) g(s) 

where the observer polynomial q(s) is an arbitrary AT th order 
stable monic polynomial, h(s) and k(s) are the unknown control 
polynomials of order N-l , and g p is the input gain [Ref 4:p. 
5]. Multiplying both sides of equation (2-7) by q(s) and 
substituting for y(t) and u(t) from the partial state 
equations (2-5) yields the closed loop dynamics 

q(s)p(s)z(t) = h(s)r(s)z(t)+k(s)p(s)z(t) + /9 _ ft v 

g p q{s)v(t ) . 

The polynomials h(s) and k(s) are defined to satisfy the 
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following condition: 



q(s)p{s)-k(s)p(s)-h(s)r(s) = — p‘(s)r(s)g(s) . (2-9) 

r a 

Equation (2-9) can be put into the form of the Diophantine 
equation [Ref. 2:p. 291] after some simple rearrangement: 



h(s)r(s)+k(s)p(s) 



g(s) (p(s)-_p*(s)r(s) ) 
r. 



( 2 - 10 ) 



where r 1 is the coefficient of the highest order term of the 
plant numerator polynomial. If the system ( p(s ) and r(s )) 
were known, then the polynomials h(s) and k(s) could be solved 
for directly using the Sylvester matrix [Ref. 2:p. 295]. The 
MATLAB subroutine FIND_HK.M was written to solve the 
Diophantine equation and return a solution for h(s) , k(s) , and 
g p . This subroutine can also be used to determine initial 
estimates of the coefficients of h(s) , k(s) , and g p for a 
linearized model of a nonlinear system. The subroutine is 
included in Appendix A. 

Assuming the estimates of h(s) and k(s) converge to the 
solution of the Diophantine equation, then the closed loop 
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partial state equations of the controlled system become: 




(2-lla) 



(2-llb) 



Eliminating the partial state variable, z(t), setting g p to 
1/r j, and dividing both sides of equation (2-lla) by q(s) 
yields the desired closed loop dynamics: 



This transfer function is identical to equation (2-6) and thus 
the closed loop system now responds to the input v(t) as would 
the reference model given by p'(s). The challenge remains to 
find the polynomials h(s) and k(s) and the gain g p that 
satisfy equation (2-9) given only the input and output data of 
a system. 

C. LIMITATIONS OF DIRECT ADAPTIVE CONTROL 

The preceding formulation of a direct adaptive control 
scheme required several assumptions which are summarized in 
Table 2-1. Of note is that the unknown plant need not be 
stable, but must be minimum phase. Because r(s) , the unknown 
plant numerator polynomial, is essentially canceled out by the 
denominator of the controller, any unstable roots of r(s) make 
the closed loop system internally unstable [Ref. 2:p. 440]. 
It should also be noted that a good estimate of the number of 



y(t) = 



i 



v(t) 



( 2 - 12 ) 



P*(s) 
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poles N and zeros M of the unknown system must be known in 
order to choose properly sized q(s) and p'(s) polynomials. 



Polynomial 


Order 


Form 


Must be stable? 


P(S) 


N 


s N +p 1 s N ‘ 1 +...p N 


No 


r(s) 


M 


r i sM1+ * • * r M 


Yes 


P*(s) 


N-M 


s N ' H +P*iS N ' M ' 1 +...pv M 


Yes 


q(s) 


N 


s N +q 1 s N ’ 1 +. . .q N 


Yes 


h(s) 


N-l 


h 1 s N * 1 +...h N 


No 


k(s) 


N-l 


k 1 s N ' 1 +...k N 


No 



Table 2-1 - Direct Adaptive Control Polynomials 



D. IMPLEMENTATION OF DIRECT ADAPTIVE CONTROL 

The traditional method of implementing the above direct 
adaptive control scheme is to use a recursive least squares 
algortihm to estimate the parameters h(s) , k(s), and g p [Ref. 
2:p. 440]. To put this problem into regression form, both 
sides of equation (2-7) are multiplied by q(s) yielding: 

g(s)u(t) = h(s)y(t)+k{s)u(t)+g p q(s)v(t) . (2-13) 

Equation (2-12) can be rewritten as: 

v(t) = p‘(s)y(t) . (2-14) 

Substituting for v(t) from equation (2-14) in equation (2-13) 
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yields 



q{s)u(t) = h(s)y(t) +k(s)u(t) +gr p q(s)p’(s)y(t) . (2-15) 



Equation (2-15) can now be written in regression form 



x(t) = d T <p(t) 



(2-16) 



where 



x(t) = q(s)u(t) ; 0 = 



V 




s^yit) 






y(t) 


*1 


? 4>(t) = 


s wl u ( t ) 


K 




u(t) 






p'(s)q(s)y(t) 



(2-17) 



Since x(t) and <P(t) are known, the least squares estimate 
of 0 is expressed as [Ref. 5: p. 324] 



S= a ;r R i n Xe.] 



(2-18) 



where 



J(Q) = Je' a(t ' r, [x(r) -Q T <t>(z) ] 2 dr. (2-19) 

0 
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The value for d is recursively estimated as 



0 (kT g +T g ) = 0 ( kT g ) + 



P(kT g +T g ) = P(kT g )~ 



P(kT e )<P(kT e ) [x(kT s )-<p T (kT a )6(kT g )) 
l+<p T (kT a )P(kT a )<p{kT a ) 

P(kT g )<P(kT a )<p T (kT a )P(kT a ) 

1 +<p T ( kT g ) P ( kT a ) <p ( kT g ) 



( 2 - 20 ) 



based on samples taken at a rate of 1/T a samples per second 
[Ref. 5:p. 325]. Although the RLS method converges to a 
correct estimate of 6, it requires extensive and intensive 
calculations, consisting of several matrix multiplications and 
scalar divisions at every iteration. Furthermore, this 
process would be difficult to implement in analog hardware 
because of the number and the nature of the required 
operations . 
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III. THE HOPFIELD NETWORK 



A. THE PROCESSING ELEMENT 

The processing element is the heart of any neural network 
and was conceived as a coarse analogy to the biological neuron 
[Ref. 6:p. NC-4]. A diagram of a typical processing element 
is shown in Figure 3-1. 




The processing element consists of three major parts: the 
input weights the summing junction; and the transfer 
characteristic g( •). The inputs x ± are either external 
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signals or signals from other processing elements. The neural 
network is defined by its structure and the values of the 
weights. The summing junction simply sums the weighted inputs 
as well as the weighted bias signal c i and passes the result 
u i to the transfer characteristic [Ref. 6:p. NC-5]. Figure 3- 
2 shows several samples of possible transfer characteristics. 




Figure 3-2 - Sample Transfer Characteristics 



The transfer characteristic is usually a monotonically 
increasing nonlinear function such as a sigmoid or a 
hyperbolic tangent [Ref 7: p. 40]. 

Neural networks are typically constructed by arranging 
processing elements in layers and interconnecting the layers. 
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B. THE HOPFIELD NETWORK 



The Hopfield network consists of a single layer of 
processing elements that are completely mutually 
interconnected [Ref. 8:pp. 96-99]. A generic Hopfield network 




Figure 3-3 - Basic Hopfield Net 

is shown in Figure 3-3. The Hopfield network uses a version 
of the processing element introduced in the previous section 
modified by the addition of an integrator at the output of 
each neuron. Thus the i th neuron in a continuous Hopfield 
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network of n elements evolves as: 



n 

d i = E t ij v j +c i < 3 ‘ 1 > 

j= i 

where v i =g(u i ), g( •) being the monotonically increasing 
transfer characteristic. Arranging the neuron states v± into 
vector V, the weight coefficients into matrix T, and the 
bias weights Ci into vector C, the equation for the entire 
Hopfield network becomes 



where 



U = TV+C 













- 




t l2 








c i 


*"21 


*22 • 


" *2„ 


• 

t 


c = 


C 2 




• 


•• t 

nn 






c 

n 



( 3 - 2 ) 



( 3 - 3 ) 



An energy function is defined for the Hopfield network of 
equation (3-2) as 



E(t) = -1v t TV-C t V ( 3 - 4 ) 

2 

Hopfield has shown that if T is negative definite and 
symmetric, then the energy function, E(t), tends to a minimum 
[Ref. 8: p. 99]. This is shown by taking the time derivative 
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of equation (3-4) along the trajectory of the network: 

£?(t) = -Lv TV- Lv r TV -C T V . (3“5) 

2 2 

Since T is symmetric, equation (3-5) is simplified to: 

E(t) = -V t TV-C t V 

= ~(V t T+C t )V (3-6) 

= -V T ( TV +C ) . 

The definition of the Hopfield network given in equation (3-2) 
is substituted in equation (3-6) yielding the time rate of 
change of the Hopfield network energy as 

f? ( t ) =-V* 17. (3“7) 

Rewriting equation (3-7) in terms of the summation definition 
of equation (3-1) yields 

£(t) = (3-8) 

i - 1 

Since g(») is monotonically increasing, its derivative, g' (•), 
must always be nonnegative. The second term in the summation, 
ii/ f is also always nonnegative. Since the derivative of Eft) 
is always nonpositive and Eft) is lower bounded, it must tend 
to a minimum. 

C. THE HOPFIELD NETWORK AS A PARAMETER ESTIMATOR 

The energy or cost function of the Hopfield network has 
been shown to have a finite minimum. Thus for any given set 
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of connection weights T and biases C the output of the 
Hopfield network V converges to a minimizing solution. This 
behavior is analogous to the recursive least squares algorithm 
whose cost function was defined in equation (2-19). The cost 
function J(6) of equation (2-19) is expanded to 



t 



i7(0) = Je" a(t ' r, x( z ) 2 dz + 

° ( 3 - 9 ) 



t 




t 


[e~ ait ~ T) <j> (z) <j> T (z) dz 
J 


0-2 


r e -a<t-t) x ( z ) 0T( r ) ( j z 

J 


0 




0 



As the first term is not a function of 6 it has no effect on 
the minimization of B and may be discarded. Thus an 
equivalent cost function for RLS estimation is: 



J(0) = 0 



t 




t 


t re' a<t " t) <^(r)<^ T (r)dr 
J 


0-2 


[e- a(t - v >x(z)<p T (z)dz 

J 


0 




0 



( 3 - 10 ) 

0 . 



A comparison between the Hopfield network energy function, 
equation (3-4), to the simplified RLS energy function, 
equation (3-10), enables certain equivalencies to be made. 
Setting the Hopfield network output V equal to 6, the 
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following equations must be true: 



T = -Je' a,t ' t, 0(r)0 r (r)dr 

t ° (3-H) 

C = fe' a(t ' t) x( r ) 0 ( r ) dr . 



Thus if the Hopfield network weight matrix and bias vector are 
set according to equation (3-11), the Hopfield network output 
will converge to the optimal estimate for 6 as would the RLS 
algorithm. 



D. THE HOPFIELD NETWORK FOR DIRECT ADAPTIVE CONTROL 

As shown above, the Hopfield network will provide an 
optimal estimate of 6 provided that T and C are formed as per 
equation (3-11). Equation (2-17) specifies the formation of 
<t>(t) as a vector of y(t), N-l derivatives of y(t), u(t), N-l 
derivatives of u(t), and the scalar p'(s)q(s)y(t). However, 
these derivatives may not be directly available to the control 
parameter identifier. Furthermore, analog differentiation is 
not a reliable operation in a real world environment as 
differentiators are highly subject to noise [Ref. 9:p. 99]. 
Rather then using the direct adaptive control equation 
directly, both sides of equation (2-15) may be operated on by 
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1/p' (s)q(s) yielding 



<?(s) 

p'(s)q(s) 



u(t) 



fr(s) 

p*(s)g(s) 

*(s) 

p‘(s)q(s) 



y(t) + 

u(t)+g p y{t). 



( 3 - 12 ) 



Equation (3-12) shows that y(t) and u(t) may be filtered by 
1/p* (s)q(s) and the resultant control parameters h(s), k(s ) , 
and g p remain the same. When y(t) and u(t) are properly 
operated on by a state space filter in controllable canonical 
form, the necessary derivative states are available without 
the need for a differentiator. 



E. CONVERGENCE AND STABILITY OF THE HOPFIELD NETWORK 

As mentioned earlier, the Hopfield network based direct 
adaptive controller will converge to an optimal estimate of 0. 
The particular transfer characteristic g(') has no effect on 
the steady state value of B unless one or more of the actual 
elements of 0 exceed the bounds of the nonlinearity. It is 
the responsibility of the designer to ensure that all 
controller parameters are within the bounds of the transfer 
characteristic . 

This implementation of the Hopfield network attempts to 
operate the neurons in their linear region. Therefore, the 
terminology Hopfield network rather than neural network is 
more appropriate to this implementation. 
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1. The effect of T and C on Hopfield network convergence 
Assuming that the Hopfield network output B does not 
exceed the bounds of the nonlinearity, then the Hopfield 
network output states behave as a linear system response to a 
step input where T is analogous to the continuous time A 
matrix and C is analogous to B. The steady state value of B 
is -T" J C and the speed of the convergence is proportional to 
the eigenvalues of T. In order to speed convergence it is 
desirable to keep the eigenvalues of T as large as possible. 

There are several methods to speed the convergence of the 
Hopfield network. The first is to be aware of the operating 
conditions of the plant in its likely area of operation. The 
rate of convergence of the algorithm is determined by the 
eigenvalues of T, and therefore by the magnitude of y(t) and 
u(t). If these signals are small enough to affect the rate of 
convergence then they can be properly scaled to increase the 
eigenvalues of T. The next method is to increase both T and 
C by a scalar positive constant X. This does not change the 
steady state value of B but it does increase the eigenvalues 
of T. In a digital simulation the use of X is almost 
unrestricted, however in an analog implementation its value is 
limited by the voltage and current capacities of the 
components and the noise level. 
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2 . The Excitation of the Input Signal 

The nature of the input signal v(t) is critical to 
parameter convergence. The signal v(t) must provide 
persistency of excitation in order to guarantee convergence of 
the parameters. Persistency of excitation is related to the 
frequency content of the signal [Ref. 5:p. 423]. For example, 
a sinusoid would not be persistently exciting for a system of 
order greater than two in the sense that it does not contain 
enough information to estimate the parameters. When an input 
signal is not persistently exciting, the eigenvalues of the T 
matrix tend to become small, slowing the Hopfield network 
convergence. Based on frequency content alone, the best input 
for persistent excitation is a white noise signal. However, 
white noise signals are not well suited for systems with small 
bandwidth. Since white noise is typically zero mean, a system 
with small bandwidth filters the white noise to a very small 
zero mean signal yielding a T matrix with very small 
eigenvalues. A superior input signal for use in these systems 
is a square wave of frequency suitable to the system since the 
square wave concentrates its energy in a finite bandwidth, 
whereas a white noise signal has an evenly spread power 
spectrum. The period of the square wave should be determined 
to ensure that it is suited to both the reference model and 
the plant. 
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F. DIGITAL SIMULATION OF A HOPFIELD NETWORK 



1. The Processing Element 

The processing element is most easily represented in 
software as a single-input single-output system. The input is 
the sum of the product of the weights and the inputs to the 
neuron added to the bias term. The output of the neuron is 
simply the output passed through the transfer characteristic. 
This work will consider four different transfer 
characteristics: 1) the sigmoid, 2) the hyperbolic tangent, 3) 
the identity (a linear neuron), and 4) a saturation 
nonlinearity. It should be noted that all four of these 
transfer characteristics are monotonically increasing. The 
subroutine that applies the nonlinearity to the neuron input 
is called SIGMOID. M and is included in Appendix A. 

2. The Hopfield Network 

Having implemented the neuron, the implementation of 
the Hopfield network is a matter of implementing the nonlinear 
differential equation that describes a Hopfield network given 
in equation (3-2). The function that iterates a Hopfield 
network over one sampling period is described below; the 
corresponding MATLAB function is included in Appendix A as 
HOPFIELD. M. If the neuron is linear (case 3 above), then U=V 
and equation (3-2) becomes a simple state space linear 
differential equation and may be simulated in a single time 
step by converting the continuous model to a discrete model 
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and iterating the discrete model by one time step. The value 
returned is the Hopfield network output at the next time step. 

The nonlinear neurons are not as simple. The MATLAB 
routine ODE45 was written to solve a nonlinear differential 
equation written in state space form [Ref 10:pp. 3-137 to 3- 
139]. This routine was modified to operate directly on the 
Hopfield network nonlinear differential equation (3-2). The 
only problem with this method is its computational complexity 
with respect to the linear method. The sigmoidal and 
hyperbolic tangent Hopfield networks were simulated in this 
manner . 

The saturation nonlinearity Hopfield network was 
implemented similarly to the linear problem. The Hopfield 
network was determined for one iteration as described for the 
linear case. Following that solution, the output was passed 
through a saturation nonlinearity. Although this was not 
strictly the solution to the nonlinear differential equation, 
it was a very close approximation. If the sampling time was 
sufficiently small, the solution arrived at through this 
simplified method was indistinguishable from that arrived at 
via the more complex ODE45 solution and there was more than a 
fifty fold savings in simulation time. 
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IV. THE HOPFIELD NETWORK FOR ADAPTIVE CONTROL 



A. THE ALGORITHM 

The implementation of the Hopfield network for direct 
adaptive control is reasonably straightforward. A schematic 
of the Hopfield network controller is shown in Figure 4-1. 




Figure 4-1 - Hopfield Net Adaptive Controller 

The first step towards the eventual implementation of the 
Hopfield network controller in hardware was a high level 
software simulation. This simulation was intended to give 
insight into the behavior of a system controlled by a Hopfield 
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network as well as to discover any potential difficulties 
inherent in this control implementation. The high level 
simulations were written in MATLAB. The next few sections 
describe the general steps in the implementation of the 
Hopfield network for direct adaptive control. 

1. Determination of the System Order 

An estimate of the system order N and the number of 
system zeros M must be determined. There are a variety of 
methods available to the designer, from complex systems 
identification tools to a simple educated guess. 

An important note in the determination of the system 
order is a caution about the modeled system zeros. If the 
optimal plant model with N poles and M zeros at a given 
operating point has unstable zeros, then the closed loop 
system will not be stable at that point. In a linear system 
this implies instability, whereas a nonlinear system may 
transit to an operating region with stable zeros and oscillate 
around the unstable operating point. One possible method to 
control this plant is to use a plant model (N and M ) of lower 
order than the one determined above. For example, if a third 
order linear plant with one unstable zero (N=3 , M=1 ) were 
controlled by this direct adaptive control algorithm, it would 
be unstable. However, if the direct adaptive controller were 
designed with the assumption that the plant was third order 
with no zeros (N=3 , M=0 ) , in some cases the closed loop system 
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would be stable because the unstable zero would never be 



canceled. Of course, this reduced form cannot model the 
unknown system as well as the full model thus reducing the 
accuracy with which the controller can follow the reference 
model. The accuracy with which a plant can follow the 
reference model is largely determined by how closely the 
reduced order optimal model matches the actual system. 

2. Determination of the Reference Model and the Observer 

The choice of the reference model, p*(s), is clearly 
critical to the closed loop system performance. Since it is 
desired that the closed loop system behave as the reference 
model to the input signal, a reference model must be chosen 
that exhibits the desired response. For this work we chose as 
reference models the class of Butterworth filters. The 
bandwidth o) 0 is chosen according to the desired speed of 
response. This reference model exhibits fast rise time with 
small overshoot. 

The observer q(s) was chosen in accordance with 
traditional control theory. A fast observer (one with very 
large poles) is not necessarily good because it will follow 
the noise as well as the signal [Ref. 5:p. 260]. As a rule of 
thumb, it is a good choice for the observer polynomial to have 
four times the bandwidth of the reference model. In this 
thesis the observer was chosen as a Butterworth polynomial 
with bandwidth 4o) 0 . 
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3. Determination of the Weight Filter Pole 



The weight filter a/ (s+a) is associated with the 
forgetting factor used in the calculation of T and C. The 
value of a is an important factor in the speed of convergence 
of the Hopfield network. For linear systems, the Hopfield 
network converges faster when a is made smaller. However, 
experimentation has shown that if a is decreased much beyond 
the slowest root of the plant, the speed of convergence will 
remain about the same. Thus a good choice for a for the 
control of a linear or nearly linear plant is the magnitude of 
the slowest root of the plant. 

The choice of a in a nonlinear plant is somewhat more 
difficult. For many nonlinear plants the above guideline for 
linear plants for choosing a is still valid. However, if the 
plant linearization changes more rapidly than the slowest root 
of the linearization, then a should be chosen to be somewhat 
larger. This increases the 'forgetting factor' of the 
Hopfield network's weights, allowing the network to 
concentrate on the more recent inputs to the system and forget 
about the older, less valid data. 

4. Determination of the Input and Output Data Filters 
The filter for the input u(t) and output y(t) data is 
determined as shown in equation (3-12) . Two vector variables, 
<P y (t) and <p u (t) , are formed by filtering y(t) and u(t) through 
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a filter with dynamics l/p.(s)q(s) in controllable canonical 
form: 





S 2 ™y f(t ) 






* y (t) = 


sy,(t) 


; *„(t> - 


su f (t) 




y f (t) 




U f (t) 



where 

y f (t) = - y(t); u f (t) = - u(t). (4-2) 

p*(s)g(s) p‘ (s)q(s) 

5. Determination of the Control Signal Filter 

The control system output from the Hopfield network 
based direct adaptive controller is given in equation (2-13). 
Both sides may be divided by q(s) to yield the control signal 

u(t) = !lifly(t)+ll±flu(t)+g v(t) (4-3) 

g(s) g(s) 

where h(s) , k(s) , and g p are the components of the parameter 
vector B. Therefore, y(t) and u(t) must be filtered by l/q(s) 
in order to generate the control signal. This is done as 
before, with two l/q(s) state space filters in controllable 
canonical form. 

B. A LINEAR SYSTEM 

Before proceeding to nonlinear systems, some basic tests 
of the Hopfield network based direct adaptive controller were 
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made on a linear system. The linear system chosen was a third 
order pitch and depth model of the AUV. The state space form 
of this system is 



4 

y 

± 



-0.07 

1 

0 



-0.04 

0 

- 0.12 



u(t) 



(4-4) 



where q is the pitch rate, y the pitch, and z the depth of the 
AUV [Ref. l:pp 27-30]. It is desired to control the depth 
state to match a reference model's response to an external 
input v(t). Figure 4-2 shows plots of the baseline run with 
the T and C filter pole a set to one radian per second, and X 
set to one. Figure 4-2a shows the reference model's and the 
actual system's response to the input signal. Figure 4-2b is 
an expanded view of the portion of Figure 4-2a outlined by the 
box. Clearly, this is not a satisfactory control as the 
output does not follow the reference model. The problem stems 
from the fact that the Hopfield network estimation of g p drops 
to nearly zero, and because g p is the coefficient that 
amplifies the input signal v(t ) , this reduces the magnitude of 
the input to the plant. The reduction of input yields a T 
matrix with extremely small eigenvalues, slowing convergence 
of the Hopfield network to a glacial pace. 

The solution to this problem is to prevent g p from falling 
below a certain threshold. The threshold value may be 
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Figure 4-2 - Hopfield control of a Linear System (g„ not 

limited, A.= l) 



determined by using system identification routines to estimate 
a model for the unknown system. Noting that g p is equal to 
l/r lf where is the highest order numerator coefficient of 
the unknown system, the threshold value may be some fraction 
of the estimate of g p or (3 !r 1 . The value of (3 is set according 
to the confidence in the estimate of r r A (3 of one implies 
absolute confidence that the actual g p is no lower than the 
estimate. A (3 of between 0.1 and 0.5 is more realistic as it 
allows some room for error in the initial estimate of g p . 
Nonlinear systems should have (3 in this range because the 
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estimate of r ; may change over time. In practice the value of 
/3 has little impact on the convergence as long as it is 
reasonably large but does not exclude the actual value of r x 
from consideration by the Hopfield network. For the next run 
of the linear system the Hopfield network estimate of g p was 
limited to 0.2g p where g p was found with the routine FIND_HK.M. 
The rest of the system parameters remain the same as the 
previous run. Figure 4 -3a shows the reference model output 
and the actual plant output due to the input signal v(t) and 
Figure 4-3b shows an expanded view of the outlined area of 
Figure 4-3a to more clearly demonstrate the plant convergence. 
It is notable that just prior to the Hopfield network 
convergence, the AUV's depth state became very large and this 
is what caused the convergence to occur so abruptly. Figure 
4-3c shows a plot of the parameter vector 6 versus time. This 
plot shows that the parameter vector converges in 
approximately 45 seconds. 

The final simulation that was run on the linear plant was 
to demonstrate the effect of increasing X to 1,000,000. The 
g p threshold was left in place as it is still required for 
rapid convergence. Figure 4-4a shows the reference model and 
the actual model and Figure 4 -4b expands the outined area of 
Figure 4 -4a and shows that the actual output does converge to 
the output of the reference model. It is again noted that the 
AUV dropped to a large depth value before the parameter vector 
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Figure 4-3 - Hopfield Control of a Linear System (g limited, 

x«u 



converged. However, the depth excursion in this case is about 
100 times less than in the previous simulation. The lower 
left plot is that of the parameter vector and shows that the 
parameter vector converges in about 10 seconds. This is a two 
fold improvement over the previous case where X was unity. 
This is as expected because the increase in X increased the 
eigenvalues of the T matrix and thus the convergence of the 
Hopfield network. Although the use of X is practically 
unrestricted in a digital simulation, it is limited by the 



33 




Figure 4-4 - Hopfield Control of a Linear System (g„ limited, 
1 = 1 , 000 , 000 ) 

available voltage and current magnitudes in an analog circuit 
implementation . 

C. THE INVERTED PENDULUM 

The inverted pendulum was used as an initial test of the 
nonlinear direct adaptive Hopfield network controller. A 
simple model will be used with the pendulum swinging on a 
stationary axis and a control torque applied at the axis. 
Figure 4-5 shows this system. The nonlinear differential 
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Figure 4-5 - Diagram of an Inverted Pendulum 
equation that describes this system is 



. Q b u 

y=— * smy- y+ (4-5) 

i ml 2 ml 2 

where b is the frictional coefficient, m is the mass of the 
pendulum and I is the length of the weightless arm. The goal 
of the control effort will be to maintain the position of the 
pendulum in a relatively upright position ( y-0 ) . In this 

upright position the system of equation (4-5) may be 
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linearized by letting sin(y)=y which yields 



1 



y(s) _ ml 2 
u(s) 




(4-6) 



1 



This linearization shows that the system order N is two, the 
number of zeros M is zero, and that the system is unstable at 
this operating point. The file PEND.M is the driver file for 
this problem and is included in Appendix A. 

For this study the following values were used: 
m=. 1 kg; b=l kg-m 2 /s; 1 = 1 m. The elements of parameter vector 
0 were initialized to 0.1. The sampling time for these 
simulations was 0.1 s. The reference model was an (N-M) th 
order Butterworth filter of bandwidth 1 r/s. As discussed 
earlier, the observer was chosen as an N th order Butterworth 
filter of frequency 4 r/s. The pole of the filter for T and 
C was 1 r/s. The input signal was a square wave of magnitude 
0.1 and frequency 0.1 r/s. The pendulum simulation was run 
with A=16*10 10 . The large value for A was required to speed 
convergence to a reasonable amount of time. Figure 4-6a shows 
a plot of the pendulum angular position in radians. Figure 4- 
6b shows a plot of the control input u(t), Figure 4-6c shows 
the time history of the Hopfield network output B(t) , and 
Figure 4-6d shows the plot of the reference signal v(t ) , the 
reference model response to v(t) , and the actual output of the 
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Figure 4-6 - Baseline Convergence of the Inverted Pendulum 



system. The convergence of B(t) is clearly evident in this 
figure. The last plot also shows that the system output 
converges to the reference model output as desired. 

Several other Hopfield network control simulations were 
attempted on the inverted pendulum with mixed results. In 
cases where the pendulum fell upside down, the Hopfield 
network had difficulty in restoring the pendulum to an upright 
position. The problem lies in the fact that this is a model 
based implementation and that the Hopfield network must be 
able to adjust its weights as quickly as the system changes 
its linearization. 
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D. THE AUTONOMOUS UNDERWATER VEHICLE 



1. AUV Fundamentals 

The model of the AUV that was used in this study is 
that of the seven foot NPS AUV [Ref 11] written in C and 
compiled for use in MATLAB [Ref. 12:pp. 124-129], The model 
is nonlinear with six degrees of freedom, 12 states, and 5 
inputs. The 12 states are: 

u 
v 
w 

p 

g 

r 
x 

y 

z 

<p 

e 

.t 

The five inputs are: the stern and bow rudder angles, the 

stern and bow planes, and the shaft RPM. A diagram of the AUV 
is shown in Figure 4-7. 

The high level control effort envisioned for the AUV 
is to control the posture of the vehicle (the posture is made 
up of the course, the x and y position, and the depth) to 



surge ( longitudinal speed) 

sway 

heave 

roll rate 
pitch rate 
yaw rate 

X 

Y 

Z (depth) 

roll 

pitch 

yaw ( course ) 
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follow a reference posture. For simplicity, the stern and bow 
actuators were treated as one control input, with the bow 
actuator receiving the negative of the signal sent to the 
stern actuator. 

2. A Control Scheme for the AUV 

Because the AUV is a multiple-input multiple-output 
nonlinear system, the control scheme is complex. The upper 
layer of control is the path-following controller which 
receives as input the posture of a reference point and the 
AUV, and outputs AUV state reference signals for use by the 
low level controller. The low level controller receives the 
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state reference signals and determines the appropriate control 
force to apply. 

a. The Path Following Algorithm 

Although the path following algorithm is not a 
major issue in this work, it is included as a possible high 
level AUV control algorithm. The control scheme implemented 
in this model was a three dimensional extension of the two 
dimensional path following algorithm described by Kanayama et 
al [Ref 13:pp. 384-389]. The path following algorithm 
compares the AUV's posture with that of a reference posture 
and determines suitable state reference signals to maintain 
the AUV on the reference path. A pictorial description of the 
postures is shown in Figure 4-8. 

The AUV is controlled by ordering the course rate, 
depth rate, and forward speed determined by the path following 
algorithm from the reference posture and the posture error. 
The path following algorithm produces three command signals 
based on the error posture between the actual vehicle and a 
reference model: 
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Depth 

Figure 4-8 - Posture Definitions 
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where the error posture is 
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It should be noted that the third command is depth change rate 
and not the heave of the vehicle and therefore it is not a 
state of the AUV model. This signal may be calculated 
analytically as 



z = -u sin0+w cos0sin$+v cos0cos <p (4-10) 

or readily estimated as 



±(kT g ) 



z(kT a )-z(kT e -T s ) 

T 



(4-11) 



where T s is the sampling interval. The goal of the control 
system is to follow the three command signals given in 
equation (4-8) by properly adjusting the control inputs. 

When the heading error is much less than one radian 
the coefficients K x , K y , K t , and K t may be interpreted and 
determined in the following manner. 1/K X is the time constant 
for the correction of the position error along the AUV’s 
longitudinal axis. K y and K t are related coefficients that 
determine a second order transfer function of the error 
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resolution across the AUV's longitudinal axis, 
differential equation for the cross range error is 



The 




( 4 - 12 ) 



Assuming the desired transfer function is to be critically 
damped, then K y and K t are related as 



Lastly, 1/K Z is the time constant for the correction of the 
depth error. 



controlled by the path-following algorithm described 
previously. The surge of the vehicle is controlled by the RPM 
command input while the depth rate is controlled by the stern 
and bow planes and the yaw rate by bow and stern rudders. The 
system order N, the number of zeros Af, and an estimate of the 
gain rj must be known for each of the three sub-systems before 
proceeding with the Hopfield network control system. These 
values were estimated by generating an input and corresponding 
output sequence for each of the sub-systems and using a system 
identification routine to estimate a linear model for the sub- 
system. The routines used to determine the system model were 
GET_MOD.M and FIND_MOD.M included in Appendix A. 




( 4 - 13 ) 



b. A Linear Model of the AUV 



There are three system outputs that are 
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(1 ) The Course Rate Controller 

An input signal for the rudders was generated 
as a square wave of magnitude 0.4 radians and frequency 0.075 
r/s with the vehicle travelling at 2 ft/s. The input signal 
was applied to the AUV model and the course rate output was 
observed. Figure 4-9a shows a plot of the input/output data 
used to determine a linear model of the rudder to course rate 
transfer function. Figures 4-9b and c present comparisons of 
the two best models. A summary of the results of the 






Figure 4-9 - Course Rate Data and Resultant Models 



system identification are shown in Table 4-1. Table 4-1 shows 
that the mean absolute state error for the first order system 
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is the second lowest. A 



third order system with 
two zeros had the 
smallest error. Figures 
4-9b and c show a 
comparison between the 
first order model and the 

Table 4-1 - System Identification 
third order model, of Course Rate Data 

Although the third order 

model is more accurate, the first order model is satisfactory. 
This reduction of model order significantly reduces the size 
of the computational problem because the size of the T matrix 
and thus the size of the Hopfield network varies as (2N+1) 2 . 
The transfer functions of the first and third order models of 
the AUV’s course rate state were estimated as 
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(4-15) 
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(2) The Depth Rate Controller 

A square wave of amplitude 0.4 radians and 
frequency 0.025 r/s was applied to the dive planes with the 
vehicle travelling at 2 ft/s. The resultant depth rate as well 
as the input signal are shown in Figure 4-10a. 






Figure 4-10 - Depth Rate Data and Resultant Models 

The output data suggests that the system is first order. The 
input and output data used with the system identification 
routine GET_M0D.M and the results are shown in Table 4-2. A 
graphic comparison of the first and fourth order model is 
shown in Figures 4-10b and c. Again the first order model has 
the second smallest summed absolute state error. A fourth 
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order model with two zeros had 
the smallest error. Although 
the fourth order model is 
slightly more accurate in terms 
of mean absolute error, the 
first order model performance is 
almost indistinguishable from 
that of the fourth order model. 

Table 4-2 - System 

The transfer functions of these Identification of Depth 

Rate Data 

models for the AUV's depth rate 
were estimated as 
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(3) The Speed Controller 

The input signal for the speed model was a 
square wave with a maximum magnitude of 480 RPM, a minimum 
magnitude of 260 RPM, and a frequency of 0.05 r/s. The AUV 
model was driven with this signal and the resultant input and 
output signals are shown in the Figures 4-lla and b 
respectively . 






Figure 4-11 - Speed Data and Resultant Models 

The system identification routine GET_M0D.M was 
again used to generate a set of linear models for the above 
input and output data sequences. The results of the system 
identification are shown in Table 4-3. 
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Since the first order model had 



the smallest mean absolute state 
error of the models considered, 
no higher order models were used 
for comparison. The estimated 
transfer function of the first 
order model of the AUV's surge 
state is 



N 


M 


Mean Abs Error 


1 


0 


0.00110 


2 


0 


0.00147 


2 


1 


0.00110 


3 


0 


0.00445 


3 


1 


00 


3 


2 


0.00128 


4 


0 


0.00468 


4 


1 


00 


4 


2 


00 



Table 4-3 - System 

Identification of Speed 
Data 



^»pe, d ( s ) 0.0007173 /A , 0 . 

- {%- lo) 

U RPH ( S ) s+0.1566 

A comparison of this model to the actual system are shown in 
Figure 4-llc. 

c. Summary of Control 

The path follower outputs command signals for the 
control of the AUV's course rate, depth rate, and surge. Three 
first order direct adaptive control Hopfield networks were 
used to provide a suitable control signal to the actuators. 
Each of the three states is modeled as a first order system 
(N=l , M=0) and the estimated models of these systems may be 
used with the routine FIND_HK.M to determine initial estimates 
of the control polynomial parameters output by the Hopfield 
network. 



49 



3 . Limitations of the Control Scheme 



The piecewise linearity of the AUV is highly dependent 
upon its forward speed [Ref. 1: pp. 25-62]. At very low 
forward speeds (i.e. at less than 1 ft/s) the control surfaces 
have little effect on the motion of the AUV making the system 
appear uncontrollable and rendering the Hopfield network 
controller and path follower ineffective. Provisions must be 
made to prevent attempted convergence of the Hopfield network 
under these poor conditions. 

4 . The AUV Control Simulation 

For the control simulation, a reference point moves 
along a prescribed path in a three dimensional space. The 
goal of the path following control is to generate three 
reference signals to keep the AUV as close to the reference 
point as possible. These three reference signals were passed 
to the Hopfield network controller in order to maintain the 
course rate, depth rate, and forward speed at the level of the 
reference signal filtered by the appropriate reference model 
l/p'(s). Thus the Hopfield network based direct adaptive 
control of the AUV requires three reference models to be 
denoted p' c (s) for course rate control, p' d (s) for depth rate 
control, and p* 8 (s) for the surge control. The determination 
of these reference models is clearly critical to system 
performance. 
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As in any real world application the available control 
force is limited. In the case of the AUV, the control force 
for course rate and depth rate control is limited by the 
physical geometry of the control surfaces. The control 
surface input is limited to about 23° or 0.4 radians of 
rotation. The shaft RPM input is limited to range between 
110 RPM and 750 RPM. The lower limit on shaft RPM prevents 
the vehicle from slowing to speeds where the control surfaces 
lose effect. These control force limits restrict the choices 
available for the reference models. If a reference model is 
too fast, the control surface is unable to meet the model and 
the system oscillates around the ordered state. The designer 
must ensure that the reference models are reasonable for the 
system to be controlled. Initially, the three models were 
chosen as Butterworth polynomials with a cut-off frequency of 
0.2 r/s. The three observer polynomials were chosen as 
Butterworth polynomials of frequency 0.8 r/s. 

The path-following parameters are also important to 
ensure accurate reference point tracking. If the path- 
follower attempts to drive the physical system beyond its 
capability to respond, the system cannot properly follow the 
path. However, if the path follower does not drive the system 
hard enough, the path following is slow to react to errors. 
The initial choice of parameters was: K x =0.5, K y =0.01, K f =0.2, 
and K=0.5. 
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The AUV was started at coordinates X=0 , Y=0 , and Z-0 . 
The reference point was started at X re{ =20, Y ref =2 0 , and Z ref =20 
generating an initial position error of 28 ft and an initial 
depth error of 20 ft. The reference point course rate, depth 
rate, and forward speed were changed according to the schedule 
in Table 4-4: 



! time 
i (s) 


course rate 
(r/s) 


depth rate 
(ft/s) 


forward speed 
(ft/s) 




0 


0 


2 


75 


n/ 20 


0.5 


2 


85 


0 


0.5 


2 


95 


tt/ 20 


0 


2 


105 


0 


-0.5 


3 


110 


0 


0 


3 


200 


n/50 


0 


2 


250 


0 


0 


2 



Table 4-4 - Reference Point Schedule 

The simulation was run and the resultant motion of the AUV in 
the XY plane is shown in Figure 4-12a. Figure 4-12b shows 
the motion of the AUV in the depth plane. Figure 4-12c shows 
the time history of the range from the vehicle to the 
reference point. These plots show that the AUV follows the 
reference point fairly well throughout the circuit. The AUV 
corrects well for the initial range error of 28 ft. The first 
two 90° turns are made at turn rates faster than the AUV model 
can match since the AUV's turn rate with full rudder 
deflection at 2 ft/s is tt/ 25 r/s. The AUV cannot keep up with 
these tight turns and the AUV's distance from the reference 
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(a) 






XY position 






Figure 4-12 - Initial AUV Path Following Simulation 



point increases. After two tight turns the AUV resolves its 
position error during the long straight leg of the path. The 
last 180° turn is at a turn rate slow enough for the AUV to 
follow and yet there is still significant cross range error, 
inferring that the cross range error coefficients are not 
properly set. The plot of the AUV's depth reveals that the 
depth control is effective with small overshoot and reasonable 
settling time. The plot of the AUV's range from the reference 
position reveals that once the initial error has been resolved 
the AUV remains within eight feet of the reference point with 
the error increasing during reference point maneuvers. 
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Although the output of the Hopfield network was 
initialized to the predetermined coefficients of the h(s) and 
k(s) polynomials and the gain q p , the linear models are not 
necessarily good representations of the AUV at varying speeds. 
Since the AUV is a nonlinear system, the Hopfield network 
output changes to adjust to new plant linearizations. Figure 
4-13 shows plots of the parameter vectors for each of the 
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Figure 4-13 - AUV Control Parameter Vector 



three controls versus time. The plot of the parameter vectors 
shows the dynamic nature of the control parameter vectors. 
The control vectors are continually adjusted to best match the 
defined reference models. 
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Figure 4-12 showed that the AUV did not satisfactorily 
follow a slow turn with the given path-following parameters. 
In order to improve the response the cross range error 
parameters were increased to K y =0.1 and K t =0.63. The path and 
initial conditions remain the same as the previous simulation. 
A plot of the system performance is shown in Figures 4-14a, b, 
and c . 




-150 -100 -50 0 50 100 150 200 



X (ft) 





Time (s) 



Time (s) 



Figure 4-14 - Improved AUV Path Following Simulation 



The AUV's path in this simulation is clearly superior 
to that of the previous run. The AUV has completely overcome 
the initial position error within 30 seconds and maintains a 
much smaller range error during the wide turn. 
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The next group of simulations are run with a different 
form of path description. The preceding path was described by 
vehicle speeds, course rates, and depth rates. Although this 
method describes a path that is fairly easy for the vehicle to 
follow, it is relatively difficult for the path planner to 
input the path. A more typical form of path description is to 
define a set of waypoints and the time of arrival at each 
waypoint. This is a far easier method for a path planner to 
output a path description. 

A set of waypoints were defined as shown in Table 4-5. 



The AUV simulation was run 












: Time 


X 


Y 


Z 


with the same constants and 


(s) 


(ft) 
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initializations as the 
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160 


20 
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previous run. The resultant 


105 


160 


70 


20 




200 


-60 


70 


20 


path of the AUV in the XY 


250 


-60 


-80 


20 



plane is shown in Figure 4- Table 4-5 - Waypoints for AUV 

Path 

15a and it shows that the AUV 

overshoots the waypoint and accumulates significant error at 
every turn. 

One method of reducing the overshoot is to reference 
the X and Y position error from a point that travels a fixed 
distance behind the reference point. This allows the AUV to 
’look ahead’ at the reference point and anticipate maneuvers. 
Several different following-ranges were tried and it was found 
that following six feet behind the reference point produced 
the best path based on mean range error. The system was 
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XY position 





Ref Path 

AUV Path 
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Figure 4-15 - AUV Traversal of 




Time (s) 

Waypoint Path 



simulated again with the waypoints from Table 4-5 and the 
resultant path in the XY plane is shown in Figure 4-16a. 

The AUV follows this path more closely because the AUV 
is able to begin turning priore to reaching the waypoints. 
However, if it is necessary for the AUV to pass through the 
waypoints this range-following modification should not be 
used. 
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X (ft) 




Figure 4-16 - AUV Traversal of Waypoint Path, Following 

Distance Six Feet 
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V. THE HOPFIELD NETWORK AS AN ELECTRONIC CIRCUIT 



The major impetus for the use of a Hopfield network for 
direct adaptive control was the ability to implement the 
Hopfield network in analog hardware. This implementation was 
a straightforward translation of the Hopfield network based 
direct adaptive controller shown in Figure 4-1. The two 
fundamental components used were an operational amplifier (op 
amp) and a four-quadrant analog voltage multiplier. Since the 
controller was intended for use in a Very Large Scale 
Integrated (VLSI) circuit the op amps and multipliers were 
designed with complementary- symmetry metal-oxide semiconductor 
(CMOS) field-effect transistor technology [Ref 9:p. 773]. 

A. A SIMPLIFIED FIRST ORDER SYSTEM CONTROLLER 

Since all three AUV controls were modeled as first order 
systems this was the system order simulated. Before 
proceeding, we considered a simplified first order direct 
adaptive controller. 

Given the uncertain first order system 

y(t) = -a y(t) +b u(t) (5-1) 

where a and b are unknown constants and the reference model 

y(t) = -p* y(t)+v(t) (5-2) 
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and equating them yields 



-a y(t)+b u (t ) = —p ’ y{t)+v(t). 



( 5 - 3 ) 



The solution of u(t) from equation (5-3) is 



u(t)=±y(t)+l(-p * y(t ) +v(t ) ) . 
b b 



( 5 - 4 ) 



When u(t) is determined by equation (5-4), then the output 
y(t) tracks the output of the reference model given by 
equation (5-2). Since a and b are unknown they must be 
estimated. 

The Hopfield network was used to estimate these parameters 
much as it was used to estimate the parameters for the higher 
order systems. Equation (5-1) was rearranged solving for u(t) 
and filtering both sides by an arbitrary first order stable 
monic polynomial q(s) 



It was seen that the coefficients a/b and 1/b present in 
equation (5-5) were the same as the unknown control 
coefficients in equation (5-4). Equation (5-5) is in the RLS 



l_u(t)=f — L_y(t)+i_L_y(t) . 



( 5 - 5 ) 



g(s) bg(s) b q(s) 
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estimation form of equation (2-15) where 



x(t) 



= — — u(t) ; 

<?(s) 



<P(t) = 



y(t) 

y(t) 




( 5 - 6 ) 



Thus a two neuron Hopfield network may be used to estimate the 
parameters a/b and 1/b necessary for the controller. A 
diagram of the first order Hopfield network controller is 
shown in Figure 5-1. 

A simulation of this implementation was conducted using 
the simulation tool TUTSIM. The code for this simulation 
implementing the flow diagram of Figure 5-1 is included in 
Appendix B. The plant to be controlled was 



y(t)=lllu(t) 
s + 5 



( 5 - 7 ) 



and the reference model was 



y(t)=_Lv(t). ( 5 - 8 ) 
s + 1 

The pole of the observer q(s) was -2 r/s, the pole of the 
filter for T and C was 100 r/s, and both outputs of the 
Hopfield network were initialized to 0.1. The reference input 
v(t) was a unit magnitude square wave of frequency 0.05 r/s. 
The plot of the input and output of the system as well as the 
Hopfield network output are shown in Figure 5-2. 
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Unknown plant 




Figure 5-1 - First Order Direct Adaptive Hopfield Controller 
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Figure 5-2a shows the input square wave v(t) and the 
system output y(t). The system output converges to the 
reference model output. Figure 5-2b shows the convergence of 
the parameters of the Hopfield network. The Hopfield network 
output converge to the vector [0.13333 0 .66667 ] T which is the 
actual value of 6. 

B. THE FIRST ORDER SYSTEM IN ANALOG HARDWARE 

Figure 5-1 was converted to an analog circuit by replacing 
the integrators, summers, and gain blocks with the 
corresponding operational amplifier circuits [Ref. 14:pp. 35- 
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125]. The multiplier blocks were replaced by Analog Devices 
internally trimmed precision IC multipliers, model number 
AD534 [Ref. 15:pp. 6-27 to 6-35]. The actual Hopfield 
network portion of the circuit worked as expected, converging 
to the expected value of 0 for a set value of T and C. 

However, when the closed loop system was run, the output 
of the Hopfield network tended to saturate. Once the Hopfield 
network output was saturated, the entire system saturated 
until the system was reset. A remedy for this problem is to 
scale the reference signal to maintain lower voltage levels in 
the system. Unfortunately, this would also tend to reduce the 
eigenvalues of the Hopfield network, slowing convergence. A 
SPICE simulation of the Hopfield network controller was 
designed to more closely analyze this problem. 

C. THE SPICE SIMULATION 

The eventual goal of this work is to generate a single 
integrated circuit (IC) that holds the aforementioned 
circuitry. Since a CMOS design is best suited to analog VLSI 
circuits, all components were designed using CMOS technology. 

1 . The CMOS Op Amp 

The CMOS op amp is a fairly common device. The 
example used in this work was chosen because of its simplicity 
as well as the availability of the SPICE parameters for the 
transistors [Ref 9:pp. 774-775]. General principles for the 
design of CMOS op amps are found in Reference 16. 
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2. The CMOS Four Quadrant Analog Voltage Multiplier 

The design of the multiplier is based on the square 
law characteristic of the current-voltage curve of the CMOS 
transistor in saturation [Ref. 17:pp. 531-532]. Figure 5-3 is 
a diagram of the CMOS multiplier. 




The resistors R 3 through R 6 are identical and transform 
the input voltages V 1 and V 2 into the transistor input signals 
Vj^/2, V 2 /2, (V 1 +V 2 )/2 for input into transistor M lf M 2 , and M 3 
respectively. Transistor M/s gate is grounded to provide a 
zero voltage reference signal. 

The four CMOS transistors are p-channel devices that 
operate in saturation. The source currents of these 



65 



transistors are 



I^K 



V. 

— -V*-V 



i 2 =k 



\2 

-±-V*-V 

2 



i 3 =k 



/ \2 

v+v. 

-1 i-v*-v 

« t 



I,K (-V-V c f 



(5-9) 



and 



WC u 

K= fi-f (5-10) 

2L 

where W and L are the length and width of the transistor gate, 
jj p is the mobility of holes, and C ox is the capacitance per 
unit area of the silicon dioxide gate. The output voltage 
referenced to V ref is 



V =-VVR (5-11) 

out 2 1 2 x ' 

where R is the resistance value of R 7 and R e . 

Thus the output of this device is a voltage difference 
proportional to the product of the two input voltages. A 
major problem with this implementation is that the voltage 
output requires a high gain differential amplifier on the 
output. Although this multiplier design is not completely 
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satisfactory, it was used in the SPICE simulation to gain 
insight into the difficulties of building an analog circuit 
direct adaptive Hopfield network controller. 

3. SPICE Simulation of the Hopfield Network 

The first step in the SPICE implementation was the 
simulation of the linear two element Hopfield network noted in 
Figure 5-1. The inputs to the system were the T matrix and 
the C vector and were set to 



-1 0.1 

0.1 -2 




(5-12) 



The theoretical solution of the Hopfield network was 
determined by solving for the steady state output of the step 
response of a state space system where A=T and B-C. The 
steady state output was calculated as [Ref. 18 :p. 688] 



y 



BB 



-T'*C = 



1.0302 
0.3015 * 



(5-13) 



The SPICE simulation was run with d initially set to [0.1 
0.1] T . A plot of the Hopfield network output is shown in 
Figure 5-4. The output of the Hopfield network converges to 
[0.9968 0.2865] T which is within five percent of the 

theoretical values shown in equation (5-13) demonstrating that 
a Hopfield network is easily implemented with analog 
components . 
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Convergence of Hopf leld Network Output 




Figure 5-4 - SPICE Hopfield Net Output 

A SPICE circuit was designed to implement Figure 5-1 
and is included in Appendix C. The major components of the 
circuit were written as subcircuits to improve the clarity of 
the code. The SPICE simulation showed that the multipliers 
representing the T and C matrices saturated causing the 
Hopfield network output to saturate as seen in the analog 
circuit mentioned earlier. One remedy to this problem is to 
scale the reference signal so that the internal signals do not 
saturate the analog devices at the cost of slowed convergence 
rate. 
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D. REMARKS ON THE ANALOG CIRCUIT IMPLEMENTATION 



The analog circuit implementation of this controller is 
theoretically possible, but the actual implementation in 
analog hardware proved more difficult. Since the multipliers 
and op amps saturate at relatively low levels, the system 
input signals may need to be scaled to ensure that all 
internal signals remain within saturation limits of the 
hardware. While the scaling process is straightforward, it 
tends to slow the rate of convergence of the Hopfield network. 

There are other possible solutions to this problem that 
may not slow the convergence rate which need to be 
investigated before a reliable analog implementation of the 
Hopfield network based direct adaptive controller can be 
completed. 
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VI. SUMMARY, CONCLUSIONS, AND RECOMMENDATIONS 



A. SUMMARY 

The stability and convergence of a direct adaptive 
Hopfield network controller was shown for linear minimum phase 
systems. This result was extended to include nonlinear 
systems that could be modeled as piecewise linear minimum 
phase systems. Simulation studies of a linear system, an 
inverted pendulum, and the NPS AUV were included to examine 
the capabilities and limitations of the Hopfield network 
control scheme. Work on a suitable path following algorithm 
was included as a possible implementation of the direct 
adaptive Hopfield network control scheme. 

The design of an electronic circuit to act as a Hopfield 
network was investigated. Computer simulations of a 
functional model of the circuit revealed no significant 
defects in the theory. However, an actual hardware model and 
a SPICE simulation of the controller did uncover several 
severe problems in the physical implementation of the circuit. 



B. CONCLUSIONS 

The simulations of the direct adaptive Hopfield network 
controller revealed the following: 
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• In a digital simulation, the Hopfield network approach to 
direct adaptive control behaves similarly to the recursive 
least squares approach. 

• The speed of convergence is highly dependent upon the 
magnitude and frequency content of the reference signal. 

• The Hopfield network controller is a suitable controller 
for use with a complex multiple input multiple output 
nonlinear system. This controller also works well within 
the framework of a higher level controller such as the 
path following algorithm. 

• The analog circuit implementation of the direct adaptive 
Hopfield network controller is feasible but is subject to 
the effects of the non-ideal analog components. 



C. RECOMMENDATIONS 

This thesis laid the foundation for further work in the 
use of Hopfield networks for direct adaptive control. There 
remains a great deal of ground uncovered including: 

• Improving the speed of convergence of the Hopfield 
network. 

• The effect of disturbance and measurement noise on the 
convergence of the Hopfield network. 

• Implementing the Hopfield network for control of non- 
minimum phase systems. 

• Optimizing the AUV path follower parameters with respect 
to the Hopfield network controller reference models. 

• Improving the analog circuit design of the direct adaptive 
Hopfield network controller to improve the circuit's 
resilience to non-ideal components. 

• Generating an analog VLSI design for the Hopfield network. 

• Writing code to provide for the automatic generation of a 
Hopfield network given the reference model, the observer 
polynomial, the system order, and the pole for the weight 
matrix filter. 
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APPENDIX A. MATLAB SOFTWARE 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FIND_HK.M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [ h , k , gp ] =f ind_hk ( a , b , ps tar , q ) 

% FIND_HK.M - USAGE: [h, k, gp] = f ind_hk(a,b,pstar,q) 

% This function determines the direct adaptive control polynomials 

% h(s) and k(s) and the input gain gp given the system denominator 

% and numerator polynomials a(s) and b(s) respectively, the 

% reference model polynomial pstar(s), and the observer polynomial 

% q(s). 

% 

% R. Scott Starsman 11-25-91 

b=rlz(b) ; 

bpstar=conv(b,pstar ) ; 
n=length(a)-l? 
m=length(b)-l; 

% Set up the Sylvester Matrix 
Sa=[a' ; zeros (n-1, 1) ] ; 

Sb=[ zeros (n-m, 1 ) ;b • ; zeros ( n-1, 1 ) ] ; 
for k=2:n 

Sa= [ Sa [ zeros (k-1, 1 ) ;a* ; zeros (n-k, 1 ) ] ] ; 

Sb=[Sb [zeros (n-m-l+k, 1) ;b' ;zeros(n-k, 1) ] ] ; 
end 

S=[Sa Sb]; % The Sylvester Matrix 

f=conv(q, a-bpstar/b( 1) ) * ; 

% Determine solution to Diophantine Equation 
hk=inv(S)*f (2:length(f ) ) ? 

% Determine controller polynomials and gain 
h=hk(n+l :2*n) ; 
k=hk( 1 :n) ; 
gp=l/b( 1 ) ; 



% Remove leading zeros 

% System Order 
% Number of Zeros 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SIGMOID. M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [ Y ] =sigmoid ( X, sigtype ) 

%SIGMOID This function subjects the input vector X to the transfer 
% characteristic designated by sigtype. The function is called 

% by: Y=sigmoid(X, sigtype) . 

% 

% sigtype specifies the transfer characteristic: 

% 1 - Sigmoid 

% 2 - tanH 

% 3 - Linear 

% 4 - Saturation 

% 

% R. S. Starsman 5-1-91 

% Copyright (c) R. S. Starsman, 1991 

% All Rights Reserved 

[n,m]=size (X) ; 
if (sigtype==l) 

Y=ones(n,m) . / (ones ( n, m) +exp( -X) ) ; 
elseif (sigtype==2 ) 

Y=(ones(n,m)-exp(-2*X) ) . / (ones (n,m) +exp ( -2*X) ) ; 
elseif ( sigtype==3 ) 

Y=X; 

elseif ( sigtype==4 ) 

Y=hard_lim(X, -1 , 1 ) ; 
end 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HOPFIELD.M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [ VI ] = hopf ield< VO , T, C , Tf , sigtype,maxval, lambda) 

%HOPFIELD Iterates a Hopfield net for Tf seconds given the initial state 
% of the Hopfield net VO, the weight matrix T, the bias vector C, 

% the sigmoid type (sigtype), the scale factor, lambda (optional), 

% and the neuron saturation level, maxval (optional). 

% 

% The matrix T must be negative definite. 

% 

% sigtype specifies the transfer characteristic: 

% 1 - Sigmoid 

% 2 - tanH 

% 3 - Linear 

% 4 - Hard lijnited saturation 

% 

% maxval is an optional parameter that specifies the upper and lower 

% bound on bounded neural outputs (sigtype=l, 2, or 4). The 

% default value is 1. 

% 

% lambda is a scale function that multiplies the T and C matrices. 

% The larger lambda is, the faster the convergence. The default 

% value is 1. 

% 

% USAGE : 

% V=hopfield(U,T,C,Tf, sigtype, maxval, lambda) 

% Set default values if necessary 
if nargin<7 
lambda=l ; 
end 

if nargin<6 
maxval=l ; 
end 

% Scale T and C by lambda 
T=lambda*T; 

C=lambda*C; 

if sigtype==3 

[ Phi , Del ] =c2d (T, C, Tf ) ; 

Vl=Phi*VO+Del; 
else 

if sigtype==4 

[Phi, Del]=c2d(T, C,Tf ) ; 

Vl=hard_lim(Phi*VO+Del, -maxval, maxval) ; % Iterate and limit output 

else 



% Linear Neuron 
% Discretize net 
% Iterate one step 



% Saturating Neuron 
% Discretize net 



% Sigmoidal or tanH Neurons 

% Routine based upon a modified version of MATLAB function ODE45.M 



% The 
alpha 
beta 



gamma 
pow = 



Fehlberg coefficients: 

- [1/4 3/8 12/13 1 1/2] ' ; 

= [ [ 1 0 0 0 

[ 3 9 0 0 

[ 1932 -7200 7296 0 

[ 8341 -32832 29440 -845 

[-6080 41040 -28352 9295 

= [ [902880 0 3953664 3855735 

[ -2090 0 22528 21970 



0 0]/4 

0 0]/32 

0 0 ] /2197 

0 0 ] /4104 

-5643 0 ] /20520 ]•; 

-1371249 277020 ]/76 18050 

-15048 -27360 J/752400 ] 



1/5; 
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trace = 0; 
tol = l.e-6; 

% Initialization 
t0=0; 

tf inal=Tf ; 
y 0=V0 ; 
t = tO; 

hmax = (tfinal - t)/5; 

hmin = (tfinal - t)/20000; 

h = (tfinal - t)/100; 

y = y0(:); 

f = y*zeros ( 1, 6 ) ; 

tout = t; 

yout = y . • ; 

tau = tol * max ( norm (y, ’inf*), 1); 

% The main loop 

while (t < tfinal) & (h >= hmin) 

if t + h > tfinal, h = tfinal - t; end 

% Compute the slopes 

% Call to neuron function sigmoid substituted here for 
% nonlinear function 

f ( : , l)=T*maxval*sigmoid(y/maxval,sigtype)+C; 
for j = 1:5 

f ( : , j+1 ) =T*maxval* sigmoid ( (y+h*f *beta ( : , j ) ) /maxval, sigtype )+C 

end 

% Estimate the error and the acceptable error 
delta = norm(h*f *gamma( : , 2 ) , * inf # ) ; 
tau = tol*max ( norm ( y, # inf * ), 1.0) ; 

% Update the solution only if the error is acceptable 
if delta <= tau 
t = t + h; 

y = y + h*f*gamma( : , 1) ; 
tout = [tout; t]; 
yout = [yout; y.*]; 

end 

% Update the step size 
if delta "= 0,0 

h = min(hmax, 0 . 8*h* (tau/delta) *pow) ; 

end 

end; 

if (t < tfinal) 

disp( * SINGUIiARITY LIKELY.*) 
t 

end 

Vl=maxval*sigmoid(yout ( length (tout ) , : ) /maxval , sigtype ) • 

end 

end 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HOPINIT «M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% HOPINIT 

% This is the initialization script file for the direct 

% adaptive control Hopfield net problem files. 

% Variable set in this file include: 

% 

% Ts - The sampling time 

% Tf - The final time 

% ZOH - the number of time steps in the zero order hold 

% sigtype - The neuron transfer characteristic (see SIGMOID. M) 

% lambda - The Hopfield net weight matrix gain 

% alpha - The pole of the filter for the T and C matrices in r/s 

% maxval - The maximum value the neurons are allowed to attain. 

% This is used for the saturating neurons 

% a - System denominator discrete time polynomial 

% b - System numerator discrete time polynomial 

% v - Reference signal 

% 

% R. Scott Starsman 11-25-91 



Ts= . 1; 

Tf =100 ; 
t=0 :Ts :Tf ; 
its=length(t ) ; 
ZOH=10 ; 
sigtype=4; 
lambda=l; 
alpha=l; 
maxval=inf ; 
rand( • normal * ) ? 



% Sampling time 
% Final time 

% # of iterations 

% Length in Ts of ZOH 

% Sigmoid type 

% T and C matrix gain 

% Filter pole for T and C matrices 



% Set up the test system (the system to be controlled is defined here) 



% The depth model of the AUV 
al= . 07 ; 
bl= . 04 ; 

v=3; % Forward speed 

k-1; 

Zd=l; % Initial depth error 



% State space model of the system 
A=[-al -bl 0;1 0 0;0 -v*bl 0]; 



B= [ 1 ; 0 ; 0 ] ; 

C= [ 0 0 1 ] ? 

D=0 ; 

x0= [ 0; 0 ; Zd] ; % 

[bc,ac]=ss2tf (A,B,C,D, 1) ; % 

[Phi, Del]=c2d(a, b,Ts ) ; % 

[b,a]=ss2tf (Phi,Del,c,d, 1) ; % 



Initial conditions 
Continuous time TF model 
Discrete time SS model 
Discrete time TF model 



% Condition a and b 
a=f liplr (a(2 : length (a) ) ) ; 
b=f liplr ( rlz (b) ) ; 



% Set system order and number of zeros 
n=length(a); % System order 

mc=length(rlz(bc) )-l; % Number of continuous time zeros 

md=length(b) -1; % Number of discrete time zeros 

% Set the number of controller parameters 
numparms=2*n+l ; 
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% Set measurement noise strength 



% Initialize system 
y=Zd*ones ( its , 1 ) ; 
u=[ zeros (2*n-l, 1 ) ; 0 . 1 ] ; 
sigmay=0; 

noise=sigmay*rand(t • ) ; 

T=zeros ( numparms , numparms ) ; 
C=zeros ( numparms , 1 ) ; 
clear theta 

theta (5,1) = . l*ones ( numparms , 1 ) ; 

% Determine the reference signal 
% 

% unit step 
%v=ones ( t ) ? 

% 

% sin wave 
%v=sin (t ) ; 

% 

% Zero 
%v=zeros (t ) ; 

% 

% Whit Noise 
%v=rand( t ) ; 

% 

% Square wave 
v=sign(sin( .2*t) ) ; 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%% H0PPR0B4.M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% H0PPR0B4.M 

% This script file simulates the direct adaptive 

% Hopfield net control of a system. The initialization 

% routine HOPINIT must be run prior to running this file. 

% 

% R. Scott Starsman 11-18-91 



rand( • normal * ) 

% Initialize copies of T and C for 

Tl=T; Cl=C ; 

tl=clock; 

ps t ar=butterw ( n-mc , 1 ) ; 
q=butterw(n,2) ? 
pstarq=conv(pstar , q) ; 
[hl,kl,gl]=find_hk(ac, be, pstar, q) ; 

% Determine th model’s response to 
model=l sim ( 1 , pstar , v ( 1 : length ( t ) ) *] 



in the ZOH 
% Set a timer 

% Determine the reference model 
% Determine the observer 

% Find the actual controller 

the input signal 

star ( length (pstar ) ) , t ) ; 



% Filter for derivatives of y and u 
[Af , Bf , Cf , Df ]=tf2ss ( 1, pstarq) ; 

Cf =eye (Af ) ; 

Df=zeros (Bf ) ; 

[Phif ,Delf ]=c2d(Af ,Bf ,Ts) ; 

ybar=zeros (2*n-mc, 1 ) ; % Initialize the filters 

ubar=zeros (2*n-mc, 1 ) ; 



% Filter for Calculating u 
Aqs=[-q(2 :n+l) ;eye(n-l) zeros (n-1, 1) ] ; 

Aq=blokdiag ( Aqs ) ; 

Bqs= [ 1 ; zeros (n-1 , 1 ) ] ; 

Bq=blokdiag(Bqs ) ; 

[Phiq, Delq]=c2d(Aq, Bq, Ts ) ; 

xq=zeros (2*n, 1) ; % Initialize the filter 

% Filter for the T and C matrices 
pole=alpha; 

A=-pole; 

B=pole; 

[PhiTC, DelTC] =c2d( A, B, Ts ) ; 

% This loop is to initialize the filters and the system 
for k=n+l:2*n 

ybar=Phif *ybar+Delf *y (k-1 ) ; 
ubar=Phif*ubar+Delf*u(k-l) ; 
xq=Phiq*xq+Delq* [y (k-1 ) ; u (k-1 ) ] ; 
u ( k ) =the ta ( : , 1 ) ' ★ [ xq ; v ( k ) ] ; 
y (k)=-a*y (k-n:k-l )+b*u(k-md-l :k-l ) ; 
end 



% Simulate the Hopfield net direct adaptive controller 
for k=2*n+l:its 

y (k)=-a*y (k-n :k-l ) +b*u (k-md-1 :k-l) ; % Iterate the system 

% Prepare phi(t) and s(t) 

ybar=Phif *ybar+Delf *y ( k-1 ) ; % Update the filtered versions 

ubar=Phif *ubar+Delf *u(k-l) ; % of y(t) and u(t) and derivs 

phi=[ybar (n-mc+1 :2*n-mc) ;ubar (n-mc+1 :2*n-mc ) ;y(k) ] ; % Form phi(t) 

s=q*ubar (n-mc :2*n-mc ) ; % Form s(t) 
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% Determine T and C and apply ZOH 

Tl=hard_lim(-DelTC*phi*phi*+PhiTC*Tl,-maxval,maxval) ; % Update T 
Cl=hard_lim(DelTC*phi*s+PhiTC*Cl,-maxval,maxval) ; % Update C 

if (k<20 ) | (rem(k,ZOH)“0) T=Tl? C=Cl; end % Apply ZOH 

% Iterate the Hopfield net 

theta( s ,k-2*n+l)=hopfield(theta( : ,k-2*n ) , T, C, Ts , lambda ,maxval, si gtype ) ; 

% Limit the estimate of gp from falling below gl/5 
if theta(numparms,k-2*n+l)>gl/5 theta(numparms,k-2*n+l)=gl/5; end 

% Filter y(t) and u(t) for calculating the control signal 
xq=Phiq*xq+Delq* [ y ( k-1 ) ; u(k-l) ] ; 

u(k)=theta( : ,k-2*n+l) * * [xq; v(k) *pstar (length (pstar) ) ] ; 
end 

etime (clock, tl) % Stop the timer 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PEND .M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% PEND .M 

% This is the MATLAB program for the Hopfield net control of 

% an inverted pendulum. 

% 

% R. Scott Starsman 11-13-91 



Ts= . 02 ; % 

Tf =100 ; % 

t=0:Ts :Tf ; 

its=length(t) ; % 

ZOH=l ; % 

ZOHstart=500; % 

sigtype=4; % 

lambda=160 000 000000; % 

pole=l; % 

maxval=inf ; 
rand( ‘normal' ) ; 

b=l ; % 

m= . 1 ; % 

1=1? % 

g=9 . 8 ; % 

n=2 ; % 

mc=0; % 

numparms=2*n+l; % 

% System initial conditions 

y= . 01*ones (its , 1 ) ; % 

yd=zeros(y); % 



u=[ zeros (2*n-l , 1 ) ? 0. 1 ] ; 

sigmay=0; 

sigmau=0 ; 

noise=sigmay*rand( t ' ) ? 

T=zeros ( numparms , nunparms ) ; 

C=zeros ( numparms , 1 ) ? 

theta ( : , 1 )=. l*ones ( numparms , 1 ) ; 



Sampling time 
Finish time 

Number of iterations 

Set the ZOH length 

Set the iteration # to start ZOH 

Transfer characteristic type 

Gain for T and C 

Pole of filter for T and C 



Pendulum damping 
Pendulum mass 
Pendulum length 
gravity 

Estimated system order 
Estimated continuous system zeros 

Number of parameters to estimate 



y 

y dot 



% This is the reference signal - square wave 
v=. 1* sign (sin ( . 5*t) )+.0*rand(t) ; 



% Hopfield Problem for upside down pendulum 



tl=clock; % 
modfreq=l; % 
pstar=butterw(n-mc , modf req) ; % 
q=butterw(n, 4*modfreq) ? % 
pstarq=conv(pstar,q) ; 



Set a timer 

Set frequency of the reference model 
The reference model 
The observer 



% Estimate the value of h(s), k(s), and gp 
[hl,kl,gl]=find_hk( [ 1 b/m/l /s 2 -g/1] , l/m/l^pstar, q) ; 



model=lsim(l,pstar,v(l:length(t) ) ,t) ? % The reference model output 



% Filter for derivatives of y and u 
[ Af , Bf , Cf , Df ]=tf2ss( l,pstarq) ; 
Cf=pstarq( length (pstarq) )*eye(Af ) ; 
Df =zeros (Bf ) ; 

[Phif ,Delf ]=c2d(Af , Bf ,Ts) ; 
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ybar^zeros ( 2*n-mc , 1 ) ; 
ubar=zeros(2*n-mc, 1) ; 

% Filter for Calculating u 
Aqs=[-q(2 :n+l) ;eye(n-l) zeros (n-1,1) ] ; 

Aq=blokdiag(Aqs) ; 

Bqs=[ 1; zeros (n-1,1) ] ; 

Bq=blokdiag(Bqs ) ; 

Cq=q( length (q) ) *eye(Aq) ; 

[Phiq,Delq]=c2d(Aq,Bq,Ts ) ; 
xq=zeros (2*n, 1 ) ; 

% Filter for the T and C matrices 
A=-pole ; 

B=pole; 

[PhiTC,DelTC]=c2d(A,B,Ts ) ; 

% Initialize the system and filters 
for k=n+l:2*n 
% The pendulum 

yd(k)=yd(k-l ) -b/m/1 A 2*yd(k-1 ) *Ts+g/l*sin (y ( k-1 ) ) *Ts; 
y(k)=y(k-l)+yd(k-l)*Ts; 

ybar=Phif *ybar+Delf *y (k-1 ) ; 
ubar=Phif*ubar+Delf*u(k-l) ; 

phi=[ybar (n-mc+1 :2*n-mc ) ;ubar (n-mc+1 : 2 *n-mc ) ;y(k-l) ] ; 
s=q*ubar (n-mc :2*n-mc ) ; 

T=hard_lim ( -DelTC*phi*phi • +PhiTC*T, -maxval , maxval ) ; 

C=hard lim(DelTC*phi*s+PhiTC*C, -maxval, maxval ) ; 
xq=(PhIq*xq+Delq* [y (k-1 ) ; u(k-l) ] ) ; 
u (k) =theta ( : , 1 ) * * [xq; v(k) ] ; 
end 

% Simulate the system 
for k=2*n+l:its 
% The pendulum 

yd(k) =yd(k-l ) -b/m/l*2*yd(k-l ) *Ts+g/l*sin(y (k-1 ) ) *Ts+u (k-1 ) /m/l~2*Ts ; 
y (^) =y (k-i >+yd(k-i> * ts; 

% Determine the phi vector and s(t) 
ybar=Phif *ybar+Delf *y (k-1 ) ; 
ubar=Phif *ubar+Delf *u(k-l ) ; 

phi= [ybar (n-mc+1 :2*n-mc ) ;ubar (n-mc+1 : 2 *n-mc ) ;y(k-l) ] ; 
s=q*ubar (n-mc :2*n-mc ) ; 

% Update T and C 

T=hard_lim ( -DelTC*phi*phi • +PhiTC*T , -maxval , maxval ) ; 
C=hard_lim(DelTC*phi*s+PhiTC*C, -maxval, maxval) ; 

% Iterate the Hopfield net 

theta ( : , k-2*n+l )=hopf ield( theta ( : , k-2*n) ,T,C,Ts, sigtype, maxval, lambda) ; 
% Prevent the estimate of gp from falling too low 

if theta (numparms,k-2*n+l )<gl/2 theta (numparms , k-2*n+l ) =gl/2 ; end 

% Filter y(t) and u(t) for the output calculations 
xq=(Phiq*xq+Delq* [y (k-1 ) ; u (k-1 ) ) ) ; 

% Determine the system control signal 
u(k)=hard__lim(theta( :,k-2*n+l) ' * [xq;v(k) ] , -30, 30 ) ; 
end 



81 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%% GET MOD «M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% GET_MOD.M 

% This file is used to determine a set of system models for 

% the input/output data held in the vectors y and u 

% respectively. These vectors are assumed to exist upon entry 

% into this routine. The routine finds a set of continuous time 

% transfer function models for varying plant order and numerator 

% order. 

% 

% R. Scott Starsman 11-25-91 



% Clear these variables for later use 
nt= [ ] ; dt=[]? results=[]; 



% Nested loop to determine the system order and number of zeros 
for n=l:5 % First to fifth order 

for m=0:n-l % # of zeros less than number of poles 

% Find the model 

[num, den , theta ]=f ind_mod(y ,u • , n,m, Ts) ; 

y2=lsim(num,den,u, t ) ; % The model’s response to the input 

e=y-y2 ; % Generate an error vector 

nt=[nt; zeros (1,6-m) num]; % Save the model numerator 

dt=[dt; zeros (1,6-n) den]; % Save the model denominator 



% Plot the actual data vs the model 
clg 

subplot (211 ) 
plot ( t,y, t,y2 , t, e ) 
subplot (212 ) 

% Plot the convergence of the model parameters 
plot ( t ( 1 : length ( theta ) ) , theta • ) 



% Save the model order, # of zeros, and summed absolute error 
results=[results;n m sum(abs (e) ) ] ; 
end 
end 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FIND_MOD.M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [num, den, theta, P] =find_mod(y,u, n,m,Ts ) 

% FIND_MOD.M - USAGE: [num, den, theta, P]=f ind_mod(y,u,n,m,Ts ) 

% This routine returns a continuous time transfer function 

% model of order n with m zeros for the input/output data y 

% and u. The time interval between data points is given by Ts • 

% 

% y - y is a column vector of system output data 

% u - u is a column vector of the system input 

% n - n is the model system order 

% m - in is the number of modeled zeros 

% Ts - Ts is the time interval between data points 

% 

% num - num is the numerator polynomial of the modeled TF 

% den - den is the denominator polynomial of the modeled TF 

% theta - theta is the vector of parameters over time. This is 

% useful for examining the convergence of the parameters 

% p - P is the error covariance matrix. 

% 

% R. Scott Starsman 11-15-91 

numparms =n+m+ 1 ; % Number of parameters to find 

theta=zeros (numparms, 1) ; % Initialize theta 

% Estimate n-1 derivatives of y and store them all in yprime 

yprime=y; 

for k=2:n 

yprime=[ [dif f (yprime ( : , 1 ) ) /Ts; 0 ] yprime ] ; 
end 

ydot=[dif f (yprime ( : , 1 ) ) /Ts; 0] ; 

% Estimate m derivatives of u and store them in uprime 

uprime=u; 

for k=2:m+l; 

upr ime= [ [ dif f ( uprime ( : , 1 ) ) /Ts ; 0 ] uprime ] ; 
end 

% Set up a matrix of y and its derivatives and u and its derivatives 
phi= [ -yprime uprime ] ; 

% Initialize error covariance matrix 
P= 1 e 8 *eye ( numparms , numparms ) ; 

% Estimate system parameters using RLS 
for k=l :length(y) -numparms 
den=l+phi(k, : )*P*phi(k, : ) • ; 

theta ( : , k+1 )=theta( : ,k)+P*phi(k, : ) * * (ydot (k) -phi(k, : ) ♦theta ( : , k) ) /den; 
P=P-P*phi(k, : ) * *phi(k, : ) ♦P/den; 
end 

% Return the numerator and denominator of the system 
den= [ 1 theta ( 1 : n , length ( theta ) ) • ] ; 
num=theta(n+l : numparms, length (theta) ) • ; 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%% AUVPATHl .M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% AUVPATHl. M 

% This file is the driver for the direct adaptive Hopfield 

% net control of the AUV. This file uses the Kanayama 

% path planning algorithm to generate state reference 

% signals for the course rate, depth rate, and surge of the 

% AUV. These state control signals are passed to the Hopfield 

% net controller for actuation of the control surfaces. 

% Initialize the system 
dt= . 25 ; 

Tf =300 ; 
t=0 : dt :Tf ; 
its=length( t ) ; 

ZOH=l ; 
lambda=l; 
sigtype=3 ; 
maxval=inf ; 
rand( • normal • ) 

% Since three states need to be controlled, three Hopfield 
% nets are reguired. Variables dealing with the course rate 
% are appended with *c*, with the depth rate are appended *d*, 

% and with the surge are appended 1 s 1 

% Set the system order and the size of the Hopfield net 
nc=l; mc=0; numparmc=2*nc+l ; 
ns=l; ms=0; numparms=2*ns+l ; 
nd=l; md=0; numparmd=2*nd+l ; 

% Initialize T, C, the input, the output, and the auv position 
Tc= zeros ( numparmc , numparmc ) ; Cc=zeros (numparmc, 1) ; gc=-5 . 6 ; 

Ts=zeros ( numparms , numpanns ) ; Cs=zeros (numparms , 1 ) ; gs=855.7; 

Td=zeros (numparmd,numparmd) ; Cd=zeros (numparmd, 1) ; gd=3.973; 
yc=0*ones (t ) ; uc=zeros (yc) ; 
ys=2*ones (t) ; us=400*ones ( t ) ; 
yd=0*ones ( t) ; ud=zeros (yd) ; 

d0=0 ; % Initial depth 

auvx=zeros ( t ) ; auvy=zeros ( t ) ; auvd=dO*ones ( t ) ; 

% Initial AUV state 

x0= [ys ( 1 ) ; 0 ; yd( 1 ) ; 0 ; 0 ; yc ( 1 ) ; 0 ; 0 ; dO ; 0 ; 0 ; 0 ] ; 

tl=clock; % Start a timer 

% The first order system models for the three subsystems 
numc=-0 .167 ;denc=[ 1 0.5]; 
nums=0 . 00065 ;dens=[ 1 0.128]; 
numd=0 . 183 ?dend= [1 0.088]; 

% The system reference models and observers 

pstarc=butterw(nc-mc, .5) ; qc=butterw(nc , 1) ; pstarqc=conv(pstarc,qc ) ; 
pstars=butterw(ns-ms, .2) ; qs=butterw(ns, . 8) ; pstarqs=conv(pstars , qs ) ; 
pstard=butterw(nd-md, .5) ; qd=butterw(nd, 1) ; pstarqd=conv(pstard, qd) ; 

% Estimate the control parameters from the system models 
[ he , kc , gc ] =f ind_hk ( dene , numc , pstarc , qc ) ; thetac= [ he ; kc ; gc ] ; 

[ hs , ks , gs ] =f ind_hk ( dens , nums , pstar s , qs ) ; thetas= [ hs ; ks ; gs ] ; 

[ hd , kd , gd ] = f ind_hk ( dend , numd , p s t ard , qd ) ; the t ad= [ hd ; kd ; gd ] ; 

% Parameters and IC*s for path follower 
% Kanayama path follower 
Kx= . 1 ; 

Ky= . 1 ; 



% Sample time 
% Final time 

% Number of iterations 
% ZOH time in numbers of time steps 
% T and C matrix gain 
% Sigmoid type 
% Neuron output limit 
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Ktheta=sqrt ( 4*Ky ) ; 

Kd=.5; 

R0=0; 
vref=2; 
wref=0; 
targddot=0 ; 

% Initial reference position 
targxyd=[20 20 20 20 20; 



in X, Y, Z 



% Critically damped system in CRE 



% Following range in ft 
% Initial reference speed 
% Initial reference course rate 
% Initial reference depth rate 



20 20 20 20 20 ? 
20 20 20 20 20 ]? 



thetaref=0? 



% Initial reference heading 



vc=wref ★ones (yc ) ? vs=vref ★ones (ys ) ? vd=targddot*ones (yd) ? 

% Filter for derivatives of y and u 
[ Af c , Bf c , Cf c , Df c ] =tf 2ss ( 1 , pstarqc ) ? 
[Phifc,Delfc]=c2d(Afc,Bfc,dt) ; 

[Af s, Bfs,Cf s, Df s]=tf2ss ( l,pstarqs ) ? 

[Phifs , Delf s ]=c2d(Af s,Bf s,dt) ? 

[Afd, Bfd, Cfd, Dfd]=tf2ss ( 1, pstarqd) ? 

[Phifd,Delfd]=c2d(Afd, Bfd,dt) ; 

ybarc=zeros (2^nc-mc, 1) ? ubarc=zeros (2^nc-mc, 1) ? 
ybars= zeros (2 ★ns -ms, 1) ; ubars=zeros (2 ★ns -ms , 1) ? 
ybard=zeros(2^nd-md, 1) ? ubard=zeros (2^nd-md, 1) ? 

% Filter for Calculating u 

Aqs= [ -qc (2 :nc+l ) ?eye (nc-1 ) zeros (nc-1, 1 ) ] ? 

Aqc=blokdiag(Aqs) ? 

Bqs= [ 1 ? zeros (nc-1 , 1 ) ] ? 

Bqc=blokdiag(Bqs) ? 

[ Phiqc , Delqc ] =c2d ( Aqc , Bqc , dt ) ; 

Aqs= [ -qs (2 :ns+l ) ?eye (ns-1 ) zeros (ns -1, 1 ) ] ? 
Aqs=blokdiag(Aqs ) ? 

Bqs=[ 1 ; zeros (ns- 1 , 1 ) ] ; 

Bqs=blokdiag(Bqs) ; 

[Phiqs, Delqs]=c2d(Aqs, Bqs,dt) ? 

Aqs= [ -qd ( 2 : nd+ 1 ) ; eye ( nd-1 ) zeros ( nd-1 f 1 ) ] ? 
Aqd=blokdiag(Aqs) ? 

Bqs= [ 1 ; zeros (nd-1 , 1 ) ] ? 

Bqd=blokdiag(Bqs) ; 

[Phiqd,Delqd]=c2d(Aqd,Bqd, dt) ? 
xqc=zeros (2^nc, 1 ) ? 
xqs=zeros(2*ns, 1) ; 
xqd=zeros(2^nd, 1) ? 



% Filter for the T and C matrices 
pole=l ; 

A=-pole ? 

B=pole ? 

[PhiTC, DelTC]=c2d(A,B / dt ) ? 

% Determine necessary length of initialization 
n=max( [nd;nc?ns] ) ? 

% Initialize the filters and the system 
for k=n+l:2*n 

ybarc=Phifc^ybarc+Delfc^yc (k-1 ) ? 
ubarc=Phifc^ubarc+Delfc*uc (k-1 ) ? 
xqc=Phiqc^xqc+Delqc^[yc(k-l) ? uc(k-l) ] ? 
uc ( k ) =hard_lim ( thetac ( : , 1 ) • ★ [ xqc ? vc ( k ) ] , - . 4 , . 4 ) ? 
ybars=Phif s^ybars+Delf s^ys (k-1 ) ? 
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ubars=Phif s*ubars+Delf s*us ( k-1 ) ; 
xqs=Phiqs*xqs+Delqs* [ys (k-1 ) ; us (k-1 ) ] ; 
us (k)=hard_lim( thetas ( : , 1 ) • ★ [xqs;vs (k) ] , 0, 750 ) ; 
ybard=Phifd*ybard+Delfd*yd(k-l ) ; 
ubard=Phifd*ubard+Delfd*ud(k-l ) ; 
xqd=Phiqd*xqd+Delqd* [yd (k-1 ) ; ud(k-l) ] ; 
ud(k)=hard_lim(thetad( 5,1) •★[xqd;vd(k) ] ,-.4, .4); 
inputs=[uc (k) ; -uc(k); ud(k); -ud(k); us(k)]; 
x0=auv2 ( xO , inputs , dt ) ; 

yc (k)=xO ( 6 ) ; ys (k)=xO ( 1 ) ; yd(k)=(xO (9 ) -dO) /dt; 
dO=xO ( 9 ) 

auvx(k)=x0(7) ; auvy (k) =xO ( 8 ) ; auvd(k)=xO (9 ) ; 
thetaref=thetaref+wref *dt; 

targxyd ( : , k+1 ) =targxyd ( : , k ) + [ cos ( thetaref ) *vref ; sin ( thetaref ) *vref ; targd 

dot]*dt; 

end 

% Simulate the AUV 
for k=2*n+l:its 

x0=auv2 (xO, inputs , dt) ; % Iterate the AUV one step 

% Pick out the course rate, depth rate, and surge 
yc(k)=x0(6) ; ys(k)=x0(l) ; yd(k)= (xO ( 9 ) -dO ) /dt; 

d0=x0 ( 9 ) ; % Keep track of the old depth 

% Update the AUV's position 

auvx(k)=x0(7) ; auvy (k) =x0 ( 8 ) ; auvd(k) =x0 ( 9 ) ; 

% Path follower definitions 

if t(k)==70 vref=2; targxyd ( 3, k) =25 ; thetaref =pi/2 ; end 

if t(k)==95 vref=2; thetaref=pi; end 

if t (k)==105 vref=3; targxyd ( 3, k) =20; end 

if t (k)==110 vref=2; end 

if t ( k) ==2 00 vref=3; thetaref =-pi/2 ; end 

if t ( k) ==250 vref=2; thetaref=0; end 

% Update the path follower 
% Range from AUV to reference point 

R(k)=sqrt ( ( targxyd ( 1 , k) -auvx(k) ) *2+( targxyd ( 2 ,k) -auvy (k) ) *2 ) ; 

% Determine the error posture 

temp= [cos (xO ( 12 ) ) sin(x0 ( 12 ) ) ; -sin(x0 ( 12 ) ) cos (xO ( 12 ) ) ] ; 
xyfollow=targxyd( 1 :2 , k) -R0* [cos (thetaref ) ; sin (thetaref ) ] ; 
e=temp* (xyfollow-[auvx(k) ;auvy(k) ] ) ; 
thetae=thetaref-xO ( 12 ) ; 

% Determine the command signals 
vc (k)=wref+vref ★ (Ky*e (2 ) +Ktheta*sin (thetae) ) ; 
vs (k)=hard lim(vref ★cos (thetae)+Kx*e ( 1 ) ,0,5); 
vd(k)=targ3dot+Kd*(targxyd(3,k)-xO(9) ) ; 

% Update the reference points posture 
thetaref=thetaref+wref *dt; 

deltarg= [ cos ( thetaref ) *vref ; sin (thetaref ) *vref ; targddot ] *dt ; 
targxyd( : , k+1 )=targxyd( : ,k)+deltarg; 

% Course Rate controller 
% Determine phi and s 
ybarc=Phifc*ybarc+Delfc*yc (k-1 ) ; 
ubarc=Phif c*ubarc+Delf c*uc (k-1 ) ; 

phi=[ybarc (nc-mc+1 :2*nc-mc ) ;ubarc (nc-mc+1 :2*nc-mc ) ;yc (k) ] ; 
s=qc*ubarc (nc-mc :2*nc-mc ) ; 
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% Determine T and C and iterate the Hopfield net 
Tc=hard_lim ( -DelTC*phi*phi • +PhiTC*Tc , -maxval , maxval ) ; 
Cc=hard_lim(DelTC*phi*s+PhiTC*Cc, -maxval, maxval) ; 
thetac( : , k-2*n+l)=. . . 

hopfield (thetac( : , k-2*n) , Tc,Cc, dt,s igtype, maxval , lambda) ; 
if thetac(numparmc,k-2*n+l)>gc/5 thetac (numparmc , k-2*n+l ) =gc/5 ; end 

% Determine the system control signal for the course rate 
xqc=Phiqc*xqc+Delqc*[yc(k-l) ; uc(k-l) ] ; 

uc(k)=thetac( : ,k-2*n+l) ' * [xqc; vc(k) *pstarc ( length (pstarc) ) ] ; 
uc(k)=hard_lim(uc(k) , -.4, .4) ; 

% Speed controller 
% Determine phi and s 
ybars=Phif s*ybars+Delf s*ys (k-1 ) ; 
ubars=Phif s*ubars+Delf s*us (k-1 ) ; 

phi=[ybars (ns-ms+1 : 2*ns-ms ) ;ubars (ns-ms+1 :2*ns-ms ) ;ys (k) ] ; 
s=qs*ubars (ns-ms : 2*ns-ms ) ; 

% Determine T and C and iterate the Hopfield net 
Ts=hard_lim( -DelTC*phi*phi ' +PhiTC*Ts , -maxval , maxval ) ; 
Cs=hard_lim(DelTC*phi*s+PhiTC*Cs , -maxval , maxval ) ; 
thetas ( : / k-2*n+l)=. . . 

hopfie Id ( thetas ( s / k-2*n) , Ts ,Cs , dt # sigtype, maxval , lambda ) ; 
if thetas (numparms / k-2*n+l )<gs/5 thetas (numparms / k-2*n+l )=gs/5; end 

% Determine the control signal for the surge 
xqs=Phiqs*xqs+Delqs* [ys (k-1 ) ; us (k-1 ) ] ; 

us (k)=thetas ( : ,k-2*n+l) • * [xqs; vs (k) *pstars ( length(pstars ) ) ] ; 
us (k)=hard_lim(us (k) , 110,750) ; 

% Depth Controller 
% Determine phi and s 
ybard=Phifd*ybard+Delfd*yd(k-l ) ; 
ubard=Phif d*ubard+Delf d*ud ( k- 1 ) ; 

phi=[ybard(nd-md+l:2*nd-md) ;ubard(nd-md+l :2*nd-md) ;yd(k) ] ; 
s=qd*ubard(nd-md:2*nd-md) ; 

% Determine T and c and iterate the Hopfield net 
Td=hard_lim( -DelTC*phi*phi • +PhiTC*Td, -maxval, maxval) ; 
Cd=hard_lim(DelTC*phi*s+PhiTC*Cd, -maxval, maxval) ; 
thetad( : ,k-2*n+l)=. . . 

hopfield (thetad( : ,k-2*n) ,Td,Cd,dt, sigtype, maxval, lambda) ; 
if thetad ( numparmd , k-2 *n+l ) <gd/5 thetad ( numparmd , k-2 *n+l ) =gd/ 5 ; end 

% Determine the control signal for the depth rate 
xqd=Phiqd*xqd+Delqd* [ yd ( k-1 ) ; ud ( k-1 ) ] ? 

ud(k)=thetad( : , k-2*n+l) ' ★ [xqd; vd ( k) *pstard ( length (pstard ) ) ] ; 
ud(k)=hard_lim(ud(k) ,-.4, .4) ; 

% Form the input vector for the AUV model 
inputs= [uc (k) ; -uc(k); ud(k); -ud(k); us(k)]; 
end 

etime (clock, tl) % stop the timer 
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APPENDIX B. TUTSIM CODE 



PROFESSIONAL VERSION OF TUTSIM 

Model File: hopfield 
Date: 11/ 23 / 1991 
Time : 22: 0 



Timing: 0.0100000 


9 


DELTA ; ,100.0000, • RANGE 


PlotBlocks and Scales 


: 






Format : 














BlockNo 


, plot-MINimum, Plot- 


MAXimum; Comment 


0,0000 , 100. 


0000 


; Time 






Yls 


26 


9 


-1# 


5000 


, 1.5000 ; y 


Y2: 


27 


9 


-1. 


5000 


, 1.5000 ; v 


Y3: 


20 


9 




5000 


, 1.5000 ; thetal 


Y4: 


21 


9 


-1. 


5000 


, 1.5000 ; theta2 


2.0000 


1 




GAI 


2 






0.0000 


2 




INT 


3 


9 


y HAT 




3 




SUM 


-1 


3 0 


y HAT dot 




4 




MUL 


3 


3 






5 




MUL 


2 


3 






6 




MUL 


3 


25 






7 




MUL 


2 


2 






8 




MUL 


2 


25 




1.0000 


9 




FIO 


-4 


9 


Til 


0.0100000 














-0.1000000 














1.0000 


10 




FIO 


-5 


. 

9 


T12 or T21 


0.0100000 














-0 1000000 














1.0000 


11 




FIO 


-7 


9 


T22 


0.0100000 














-0.1000000 














1.0000 


12 




FIO 


8 


9 


C2 


0.0100000 














0.1000000 














1.0000 


13 




FIO 


6 


9 


Cl 


0.0100000 














0.1000000 
















14 




MUL 


20 


9 






15 




MUL 


20 


10 






16 




MUL 


21 


10 






17 




MUL 


21 


11 






18 




SUM 


13 


14 16 






19 




SUM 


12 


15 17 




0.1000000 


20 




INT 


18 


9 


thetal 


0.1000000 


21 




INT 


19 


9 


theta2 




22 




MUL 


20 


29 






23 




MUL 


21 


26 






24 




SUM 


22 


23 


u 


1.0000 


25 




FIO 


24 


9 


U HAT 


0.5000000 














0.0000 
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1.5000 

0.2000000 

0.0000 


26 


FIO 


28 


? y 




27 


SGN 


31 


? v 


1.0000 


28 


GAI 


24 






29 


SUM 


-26 27 






30 


GAI 


26 




0.0500000 

1.0000 


31 


FRQ 
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APPENDIX C. SPICE SOFTWARE 



FIRST ORDER HOPFIELD NET CONTROLLER 
* IMPORTANT NODES: 



★ 


1 


-> 


y 






★ 


2 


-> 


yHAT 






★ 


3 


-> 


yHATDOT 






★ 


i 

i 


9 


io ! -> ! 


Til 


T12 


★ 


i 

i 


10 


ii ! -> i 


T21 


T22 


★ 


i 

i 


12 


! -> ! ci 


i 

i 




★ 


i 

i 


13 


! -> j c2 


i 





* 20 -> THETAl 

* 21 -> THETA2 

* 24 -> U 

* 25 -> UHAT 
» 27 -> v 

XI 24 101 102 1 SYSTEM 
X2 1 101 102 2 3 YFILTER 
X3 3 3 101 102 4 MULT 

X4 2 3 101 102 5 MULT 
X5 2 2 101 102 6 MULT 
X6 3 25 101 102 7 MULT 
X7 2 25 101 102 8 MULT 
X8 4 101 102 9 TFILTER 
X9 5 101 102 10 TFILTER 
X10 6 101 102 11 TFILTER 

XII 7 101 102 12 CFILTER 
X12 8 101 102 13 CFILTER 
X13 9 20 101 102 14 MULT 
X14 10 20 101 102 15 MULT 
X15 10 21 101 102 16 MULT 
X16 11 21 101 102 17 MULT 

X17 12 14 16 101 102 18 SUMMER3 

X18 13 15 17 101 102 19 SUMMER3 

X19 18 101 102 20 INTEG 

X20 19 101 102 21 INTEG 

X21 20 26 101 102 22 MULT 

X22 1 21 101 102 23 MULT 

X23 22 23 101 102 24 SUMMER2 

X24 28 101 102 24 INVERTER 

X25 24 101 102 25 UFILTER 

X26 27 1 101 102 26 DIFF 

W 27 0 PULSE (-5 5 .1 0 0 4 8) 

VDD 101 0 DC 12V 

VSS 102 0 DC -12V 

.TRAN 100MS 20S 

.PLOT TRAN V(l) V(27) 
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. SUBCKT SYSTEM 3124 

* u Vdd Vss y 

* FIRST ORDER SYSTEM — > y/u=G/(Ts+l) 

* G= (R1+R2 ) /Rl ; T=R3*C1 

Rl 6 0 100K 

R2 4 6 5 OK 
R3 3 5 200K 
Cl 5 0 10 

XI 5 6 1 2 4 MOPAMP 
•ENDS SYSTEM 

•SUBCKT CFILTER 3124 

* cin Vdd Vss cout 

* FILTER FOR C VECTOR — > COUt/Cin=G/ (Ts+1 ) 

* G= 1 ; T=R3*Cl 
Rl 4 0 100K 

R3 3 5 100K 
Cl 5 0 1U 

XI 5 4 1 2 4 MOPAMP 
•ENDS CFILTER 

•SUBCKT TFILTER 3124 

* Tin Vdd Vss Tout 

* FILTER FOR T MATRIX — > Tout /Tin=-G/ ( Ts+1 ) 

* G= 1 ; T=R1*C1; Rl=R2=R3 ; R4=2*Rl 

Rl 3 7 100K 

R2 6 7 100K 

R3 4 6 200K 

R4 5 0 200K 

Cl 7 0 1U 

C2 4 6 .5U 

XI 5 6 1 2 4 MOPAMP 

•ENDS TFILTER 

•SUBCKT YFILTER 11 1 2 8 3 

* y Vdd Vss yHAT yHATDOT 

* FILTER TO DERIVE yHAT AND yHATDOT — > yHAT/y=l/ (Ts+1 ) ; 
y HATDOT/y=S / ( TS+ 1 ) 

* T=R2/R1=R4/R3 
Rl 11 12 10K 

R2 12 0 20K 

R3 9 8 10K 

R4 9 3 2 OK 

XI 12 9 1 2 3 MOPAMP 

X2 3 1 2 5 INTEG 

X3 5 1 2 8 INVERTER 

•ENDS YFILTER 

•SUBCKT UFILTER 3124 

* Uin Vdd Vss Uout 

* FILTER FOR C VECTOR — > Uout/Uin=G/ (Ts+1 ) 

* G=l; T=R3*Cl 
Rl 4 0 500K 

R3 3 5 50 OK 
Cl 5 0 1U 

XI 5 4 1 2 4 MOPAMP 
•ENDS UFILTER 
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4 

vout 



• SUBCKT INVERTER 312 

♦ Vin Vdd Vss 

★ INVERTER — > Vout=-G*Vin 

G=R2/R1; R3=R1 | ] R2 



R1 

R2 

R3 

XI 



3 6 100K 

4 6 100K 

5 0 5 OK 

5 6 12 4 



MOPAMP 



.ENDS INVERTER 

.SUBCKT DIFF 3 4 

* X Y 

* DIFFERENCE CIRCUIT 
R1 3 6 100K 

100K 
100K 
100K 
12 5 



1 2 

Vdd Vss 



5 

X-Y 



R2 

R3 

R4 

XI 



MOPAMP 



.ENDS DIFF 



.SUBCKT SUMMER2 3 4 

★ X Y 

♦ 2 INPUT SUMMER CIRCUIT 



1 

Vdd 



2 

Vss 



5 

-<X+Y) 



R1 

R2 

R3 

R4 

XI 



100K 
100K 
100K 
33K 
12 5 



MOPAMP 



.ENDS SUMMER2 

.SUBCKT SUMMER3 3 4 

* X Y 

* 3 INPUT SUMMER CIRCUIT 
Rl 3 8 100K 



12 6 

Vdd Vss -(X+Y+Z) 



R2 

R3 

R4 

R5 

XI 



8 100K 
8 100K 
8 100K 
0 25K 
8 12 6 



MOPAMP 



.ENDS SUMMER3 

.SUBCKT INTEG 3124 

★ Vin Vdd Vss Vout 

★ INTEGRATOR — > VOUt/Vin=G/s 

★ G=l/ (Rl*Cl ) ; R2=R1 

Rl 3 5 lMEG 

R2 6 0 lMEG 
Cl 4 5 1U 

XI 6 5 1 2 4 MOPAMP 
.ENDS INTEG 



.SUBCKT MULT 3 4 1 2 12 

★ X Y Vdd Vss Vout 

★ANALOG VOLTAGE MULTIPLIER FROM HONG & MELCHIOR (ELEC LTTRS 6 JUN 85) 

Rl 3 5 5K 

R2 0 5 5K 

R3 4 6 5K 

R4 0 6 5K 

R5 4 7 5K 

R6 3 7 5K 

Ml 1 5 8 1 CMOSP L=32U W=16U 
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M2 1 6 8 1 CMOSP L=32U W=16U 

M3 1 7 9 1 CMOSP L=32U W=16U 

M4 1 0 9 1 CMOSP L=32U W=16U 

XI 9 8 1 2 10 MOPAMP 

R7 11 9 10K 

R8 10 8 10K 

El 12 0 10 11 297 

ROUT 12 0 1M 

V2MINUS 11 0 DC -10 

.ENDS MUI.T 

. SUBCKT COPAMP3 34125 



★ 








V+ V- 


Vdd Vss 


Vout 


♦CMOS 


OPAMP FROM 


MICROELECTRONIC 


CIRCUITS (SEDRA AND SMITH) 


Ml 9 


4 


6 


1 


CMOSP 


L=8U W=120U 


M2 8 


3 


6 


1 


CMOSP 


L=8U W=120U 




M3 9 


9 


2 


2 


CMOSN 


L= 10U 


W=50U 




M4 8 


9 


2 


2 


CMOSN 


L=10U 


W=50U 




M5 6 


7 


1 


1 


CMOSP 


L=10U 


W=150U 




M6 5 


8 


2 


2 


CMOSN 


L=10U 


W=100U 




M7 5 


7 


1 


1 


CMOSP 


L=10U 


W=150U 




M8 7 


7 


1 


1 


CMOSP 


L= 1 0U 


W=150U 




Cl 8 


10 


1 0PF 








R4 10 




5 


10K 








M9 7 


7 


2 


2 


CMOSN 


L=2 50U W=5U 




.ENDS 


COP AMP 3 








.SUBCKT 


MOPAMP 1 


2 


98 99 


3 


★ 








V+ V- 


Vdd Vss 


Vout 



RIN 1 2 6 MEG 
ROUT 4 3 75 



El 0 4 1 2 200K 

* R1 AND R2 ARE USED TO MAKE THE OPAMP IDEAL MODEL PIN COMPATIBLE 

* WITH THE CMOS OPAMP 
R1 98 0 10M 

R2 99 0 10M 
.ENDS MOPAMP 

•MODEL CMOSN NMOS LEVEL=2 . 00000 LD=0.6U TOX=8.5E-8 

+NSUB=1E+16 VTO= 1 JS=lE-6 UO=750 .000 UEXP=. 14 UCRIT=5E4 UTRA=0 PB=.7 
+VMAX=5E4 XJ=lU CGBO=2E-10 

+RSH=15 CGSO=4E-10 CGDO=4E-10 CJ=4E-04 MJ=2 CJSW=8E-10 MJSW=2 

•MODEL CMOSP PMOS LEVEL=2 . 00000 LD=0.6U TOX=lE-7 NSUB=2E+15 VTO=-l JS=lE-6 

+U0=100.000 UEXP= .03 UCRIT=lE4 UTRA=0 PB=.7 VMAX=3E4 XJ=.9U CGBO=2E-10 

+RSH=75 CGSO=4E-10 CGDO=4E-10 CJ=1.8E-04 MJ=2 CJSW=6E-10 MJSW=2 

.END 
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